An improved computational electromagnetics process and system

ABSTRACT

A computational electromagnetics system has a memory and at least one processor configured to execute a computational electromagnetics process. The process includes solving a pair of differential equations for only two variables, the two variables representing a pair of scalar potentials, and determining at least one of an electric field and a magnetic field from the pair of scalar potentials.

TECHNICAL FIELD

The present invention relates to an improved computationalelectromagnetics process and system, in particular to a more efficientprocess for generating numerical solutions to problems inelectromagnetics defined by partial differential equations.

BACKGROUND

Computational electromagnetics plays a paramount role in the modelingand understanding of electromagnetic phenomena, and generally involvesaccurately modeling electromagnetic problems through robust numericalmethods. Computational electromagnetics methods are used to calculateelectromagnetic fields for a wide array of different applications,including not only the design of many types of physical apparatus,including for example integrated circuits, aircraft, wireless systems,telecommunications and radar systems, but also for advanced medicalimaging, diagnosis and treatment methods such as detecting and treating(hyperthermia) tumors, for example. In many of these applications (e.g.,microwave-based imaging, radar beamforming or tomography), there is ageneral need to reduce the calculation time as much as possible.

The most popular and efficient computational electromagnetics methodsare finite difference methods (FDM), finite element methods (FEM),moment methods (MOM) and meshless methods, and state of the artcommercial software packages, such as HFSS, CST MS, ADS, COMSOL, XFdtdand OmniSim FDTD, are based on one or more of those computationalmethods.

Despite the considerable power of existing computationalelectromagnetics methods in their ability to generate accurate solutionsto electromagnetic problems, a continuing difficulty of these methods isthe substantial (real and processor) time required to generate accuratesolutions to the underlying partial differential equations. One approachto reducing the (wall clock) time required to generate solutions tocomplex computational electromagnetics problems is to use asupercomputer. For example, the Computational ElectromagneticsLaboratory of NASA provides access to a computer cluster with 476processors and 1.95 TB of random access memory dedicated tocomputational electromagnetics. Nevertheless, several days ofcomputation may be required to simulate the effect of electromagneticfields on a large scale structure such as an aircraft.

Consequently, there is a general desire or need to reduce thecomputation times as much as practically possible, and for someapplications (e.g., ‘real-time’ imaging and monitoring) there would bedistinct advantages if such calculations could be performed in timesapproaching ‘real-time’, where the delays are in substance negligible.

In view of the above, many attempts have been made to increase thetime-efficiency of the aforementioned classical methods. However,despite these ongoing efforts, even today's state of the art methodsstill require inconveniently long times to generate solutions toelectromagnetic problems, especially in the case of large scaleproblems.

It is desired, therefore, to address or alleviate one or moredifficulties of the prior art, or to at least provide a usefulalternative.

SUMMARY

In accordance with some embodiments of the present invention, there isprovided a computational electromagnetics process for determiningelectromagnetic fields, the process being executed by at least oneprocessor of a data processing system, and including the steps of:

-   -   (i) solving a pair of differential equations for only two        variables, the two variables representing a pair of scalar        potentials; and    -   (ii) determining at least one of an electric field and a        magnetic field from the pair of scalar potentials determined at        step (i).

In some embodiments, step (ii) includes at least one of:

-   -   (iii) determining the electric field {right arrow over (E)}        according to:

{right arrow over (E)}=jωψ∇φ; and

-   -   (iv) determining the magnetic field {right arrow over (B)}        according to:    -   where φ and ψ represent the scalar potentials, ω represents        angular frequency, and j=√{square root over (−1)}.

In some embodiments, the differential equations are of the form:

∇² ψ+k ²ψ=0

∇²φ=0

-   -   where k=ω√{square root over (εμ)} is the wavenumber, ε        represents permittivity, and μ represents permeability.

In some embodiments, the differential equations are of the form:

∇² ψ+k ²ψ=0

∇² v+k ² v=0

-   -   where k=ω√{square root over (ωμ)} is the wavenumber, ε        represents permittivity, μ represents permeability, and v        represents scalar electric potential, and wherein

$v = {- {\frac{j\; {\omega\phi\psi}}{2}.}}$

In some embodiments, the differential equations are of the form:

∇²ψ + (2k² − γ²)ψ = 0 ∇²ϕ − γ²ϕ = 0 or ∇²ψ + (2k² − γ²)ψ = 0${{\nabla^{2}v} + {k^{2}v}} = {- \frac{\rho + \rho_{J}}{ɛ}}$

-   -   where φ and ψ represent the scalar potentials,

${\gamma = {{- \frac{j}{2}}k}},$

k=ω√{square root over (εμ)} wavenumber, ε represents permittivity, μrepresents permeability, v represents scalar electric potential, andv=−jωφψ.

In accordance with some embodiments of the present invention, there isprovided a computational electromagnetics process executed by at leastone processor of a data processing system, the process including atleast one of:

-   -   determining an electric field {right arrow over (E)} according        to:

{right arrow over (E)}=jωψ∇φ; and

-   -   (ii) determining a magnetic field {right arrow over (B)}        according to:

{right arrow over (B)}∇φ×∇ψ;

-   -   where ω represents angular frequency, j=√{square root over        (−1)}, and φ and ψ represent scalar potentials satisfying a pair        of differential equations for only two variables, the two        variables representing the scalar potentials.

In some embodiments, the pair of differential equations are of the form:

∇² ψ+k ²ψ=0

∇²φ=0

-   -   where k=ω√{square root over (εμ)} is the wavenumber, ε        represents permittivity, and μ represents permeability.

In some embodiments, the pair of differential equations are of the form:

∇² ψ+k ²ψ=0

∇² v+k ² v=0

-   -   where k=ω√{square root over (εμ)} is the wavenumber, ε        represents permittivity, μ represents permeability, and v        represents scalar electric potential, and wherein

$v = {- {\frac{j\; {\omega\phi\psi}}{2}.}}$

In some embodiments, the pair of differential equations are of the form:

∇²ψ + (2k² − γ²)ψ = 0 ∇²ϕ − γ²ϕ = 0 or ∇²ψ + (2k² − γ²)ψ = 0${{\nabla^{2}v} + {k^{2}v}} = {- \frac{\rho + \rho_{J}}{ɛ}}$

-   -   where k=ω√{square root over (εμ)} is the wavenumber,

${\gamma = {{- \frac{j}{2}}k}},$

ε represents permittivity, μ represents permeability, v representsscalar electric potential, and v=−jωφψ.

In accordance with some embodiments of the present invention, there isprovided at least one computer-readable storage medium having storedthereon executable instructions that, when executed by at least oneprocessor of a data processing system, cause the at least one processorto execute the computational electromagnetics process of any one of theabove processes.

In accordance with some embodiments of the present invention, there isprovided at least one computer-readable storage medium having storedthereon executable instructions that, when executed by at least oneprocessor of a data processing system, cause the at least one processorto execute a computational electromagnetics process, including the stepsof:

-   -   (i) solving a pair of differential equations for only two        variables, the two variables representing a pair of scalar        potentials; and    -   (ii) determining at least one of an electric field and a        magnetic field from the pair of scalar potentials determined at        step (i).

In some embodiments, step (ii) includes at least one of:

-   -   (iii) determining the electric field E according to:

{right arrow over (E)}=jωψ∇φ; and

-   -   (iv) determining the magnetic field {right arrow over (B)}        according to:

{right arrow over (B)}=∇φ×∇ψ;

-   -   where φ and ψ represent the scalar potentials, ω represents        angular frequency, and j=√{square root over (−1)}.

In some embodiments, the differential equations are of the form:

∇² ψ+k ²ψ=0

∇²φ=0

-   -   where k=ω√{square root over (εμ)} is the wavenumber, ε        represents permittivity, and μ represents permeability.

In some embodiments, the differential equations are of the form:

∇² ψ+k ²ψ=0

∇² v+k ² v=0

-   -   where k=ω√{square root over (εμ)} is the wavenumber, ε        represents permittivity, μ represents permeability, and v        represents scalar electric potential, and wherein

$v = {- {\frac{j\; {\omega\phi\psi}}{2}.}}$

In some embodiments, the differential equations are of the form:

∇²ψ + (2k² − γ²)ψ = 0 ∇²ϕ − γ²ϕ = 0 or ∇²ψ + (2k² − γ²)ψ = 0${{\nabla^{2}v} + {k^{2}v}} = {- \frac{\rho + \rho_{J}}{ɛ}}$

-   -   where φ and ψ represent the scalar potentials,

${\gamma = {{- \frac{j}{2}}k}},$

k=ω√{square root over (εμ)} is the wavenumber, ε representspermittivity, μ represents permeability, v represents scalar electricpotential, and v=−jωφψ.

In accordance with some embodiments of the present invention, there isprovided a computational electromagnetics system having a memory and atleast one processor configured to execute a computationalelectromagnetics process, including the steps of:

-   -   (i) solving a pair of differential equations for only two        variables, the two variables representing a pair of scalar        potentials; and    -   (ii) determining at least one of an electric field and a        magnetic field from the pair of scalar potentials determined at        step (i).

In some embodiments, step (ii) includes at least one of:

-   -   (iii) determining the electric field {right arrow over (E)}        according to:

{right arrow over (E)}=jωψ∇φ; and

-   -   (iv) determining the magnetic field {right arrow over (B)}        according to:

{right arrow over (B)}=∇φ×∇ψ;

-   -   where φ and ψ represent the scalar potentials, ω represents        angular frequency, and j=√{square root over (−1)}.

In some embodiments, the differential equations are of the form:

∇² ψ+k ²ψ=0

∇²φ=0

-   -   where k=ω√{square root over (εμ)} is the wavenumber, ε        represents permittivity, and μ represents permeability.

In some embodiments, the differential equations are of the form:

∇² ψ+k ²ψ=0

∇² v+k ² v=0

-   -   where k=ω√{square root over (εμ)} is the wavenumber, ε        represents permittivity, μ represents permeability, and v        represents scalar electric potential, and wherein

$v = {- {\frac{j\; {\omega\phi\psi}}{2}.}}$

In some embodiments, the differential equations are of the form:

∇²ψ + (2k² − γ²)ψ = 0 ∇²ϕ − γ²ϕ = 0 or ∇²ψ + (2k² − γ²)ψ = 0${{\nabla^{2}v} + {k^{2}v}} = {- \frac{\rho + \rho_{J}}{ɛ}}$

-   -   where φ and ψ represent the scalar potentials,

${\gamma = {{- \frac{j}{2}}k}},$

k=ω√{square root over (εμ)} is the wavenumber, ε representspermittivity, μ represents permeability, v represents scalar electricpotential, and v=−jωφψ.

In accordance with some embodiments of the present invention, there isprovided a computational electromagnetics process for determiningelectromagnetic fields, the process being executed by at least oneprocessor of a data processing system, and including the steps of:

-   -   (i) solving a pair of differential equations for only two        variables, the two variables representing a pair of scalar        potentials; and    -   (ii) determining at least one of an electric field and a        magnetic field from the pair of scalar potentials determined at        step (i);

wherein the differential equations are forms of the Schrödinger equationof quantum physics.

The differential equations may be of the form:

∇² ψ+k ²ψ=0

∇²φ=0

where φ and ψ are scalar variables;

and wherein step (ii) includes at least one of:

-   -   (iii) determining an electric field according to:

{right arrow over (E)}=jωψ∇φ; and

-   -   (iv) determining a magnetic field according to:

{right arrow over (B)}∇φ×∇ψ;

-   -   where ω represents angular frequency, and j=√{square root over        (−1)}.

In accordance with some embodiments of the present invention, there isprovided a computational electromagnetics process executed by at leastone processor of a data processing system, the process including atleast one of:

-   -   (i) determining an electric field according to:

{right arrow over (E)}=jωψ∇φ; and

-   -   (ii) determining a magnetic field according to:

{right arrow over (B)}=∇φ×∇ψ;

-   -   where ω represents angular frequency, j=√{square root over        (−1)}, and φ and ψ are scalar variables satisfying:

∇² ψ+k ²ψ=0

∇²φ=0

In accordance with some embodiments of the present invention, there isprovided at least one computer-readable storage medium having storedthereon executable instructions that, when executed by at least oneprocessor of a data processing system, cause the at least one processorto execute a computational electromagnetics process, including at leastone of:

-   -   (i) determining an electric field according to:

{right arrow over (E)}=jωψ∇φ; and

-   -   (ii) determining a magnetic field according to:

{right arrow over (B)}=∇φ×∇ψ;

-   -   where ω represents angular frequency, j=√{square root over        (−1)}, and φ and ψ are scalar variables satisfying:

∇² ψ+k ²ψ=0

∇²φ=0

In accordance with some embodiments of the present invention, there isprovided a computational electromagnetics system having a memory and atleast one processor configured to execute a computationalelectromagnetics process, including at least one of:

-   -   (i) determining an electric field according to:

{right arrow over (E)}=jωψ∇φ; and

-   -   (ii) determining a magnetic field according to:

{right arrow over (B)}=∇φ×∇ψ;

-   -   where ω represents angular frequency, j=√{square root over        (−1)}, and φ and ψ are scalar variables satisfying:

∇² ψ+k ²ψ=0

∇²φ=0

BRIEF DESCRIPTION OF THE DRAWINGS

Some embodiments of the present invention are hereinafter described, byway of example only, with reference to the accompanying drawings,wherein:

FIG. 1 is a block diagram of a computational electromagnetics system inaccordance with an embodiment of the present invention;

FIG. 2 is a flow diagram of a computational electromagnetics processexecuted by the computational electromagnetics system;

FIG. 3(a) is a representation of a simple dipole antenna operating at afrequency of 1 GHz;

FIG. 3(b) is a corresponding set of three colour representations of thepolar distributions of electric potential (voltage) in the xy-plane(left), yz-plane (middle) and the zx-plane (right);

FIG. 3(c) is a corresponding set of three colour representations of thepolar distributions of the norm of the scalar potential yi in thexy-plane (left), yz-plane (middle) and the zx-plane (right);

FIGS. 3(d) to (g) are corresponding pairs of colour plots representingthe simulated normalized electric ((d), (e)) and magnetic ((f), (g))field radiation patterns generated by the dipole antenna of FIG. 3(a)operating at 1 GHz, as simulated using the classical (plots (d), (f))and two scalar potential (plots (e), (g)) methods in the xy-plane (wavefront, top row), the yz-plane (H-plane, middle row), and zx-plane(E-plane, bottom row);

FIG. 4(a) shows a spherical scatterer arranged in front of the dipoleantenna of FIG. 3(a);

FIG. 4(b) is a corresponding set of three colour representations of thepolar distributions of electric potential (voltage) in the xy-plane(left), yz-plane (middle) and the zx-plane (right);

FIG. 4(c) is a corresponding set of three colour representations of thepolar distributions of the norm of the scalar potential yi in thexy-plane (left), yz-plane (middle) and the zx-plane (right);

FIGS. 4(d) to (g) are corresponding pairs of colour plots representingthe simulated normalized electric ((d), (e)) and magnetic ((f), (g))field radiation patterns generated by the dipole antenna of FIG. 3(a)operating at 1 GHz, as simulated using the classical (plots (d), (f))and two scalar potential (plots (e), (g)) methods in the xy-plane (wavefront, top row), the yz-plane (H-plane, middle row), and zx-plane(E-plane, bottom row);

FIG. 5(a) shows eight antennas identical to the dipole antenna of FIG.3(a) arranged about a model of a human head (see text for details);

FIG. 5(b) is a corresponding set of three colour representations of thepolar distributions of electric potential (voltage) in the xy-plane(left), yz-plane (middle) and the zx-plane (right);

FIG. 5(c) is a corresponding set of three colour representations of thepolar distributions of the norm of the scalar potential yz in thexy-plane (left), yz-plane (middle) and the zx-plane (right);

FIGS. 5(d) to (g) are corresponding pairs of colour plots representingthe simulated normalized electric ((d), (e)) and magnetic ((f), (g))field radiation patterns generated by the dipole antenna of FIG. 3(a)operating at 1 GHz, as simulated using the classical (plots (d), (f))and two scalar potential (plots (e), (g)) methods in the xy-plane (wavefront, top row), the yz-plane (H-plane, middle row), and zx-plane(E-plane, bottom row).

DETAILED DESCRIPTION

As known by those skilled in the art, existing computationalelectromagnetics processes and systems all solve the governing partialdifferential equations of electromagnetic problems using the auxiliaryscalar potential v and vector electromagnetic potential {right arrowover (A)}. The electric and magnetic fields in general structures arethen calculated according to these four variables ({right arrow over(A)} being a vector with 3 spatial components or ‘variables’). Insource-free problems (e.g., cavities), the electric field and themagnetic field are each directly derived using the above methods bysolving a 3-variable field problem (electric and magnetic fields beingvectors with three spatial components each). Although these approachesare able to generate accurate results, the computational resourcesrequired are undesirably high.

The described embodiments of the present invention constitute asubstantial advance over the prior art by providing a fundamentally newapproach to solving computational electromagnetics (EM), requiring onlytwo variables and significantly reducing the computational resourcesrequired to generate accurate solutions to electromagnetic problems(e.g., around 70% reduction in computation time).

As known by those skilled in the art, electromagnetic theory is governedby the four Maxwell equations, together with the continuity condition asthe fifth equation. In general, those equations are usually solved forany EM problem to find the electric field ({right arrow over (E)}) andmagnetic flux density ({right arrow over (B)}) It is difficult todirectly calculate, in an analytical or numerical way, {right arrow over(E)} and {right arrow over (B)} due to their complex vectorrepresentation in Maxwell equations. For example, in the case ofisotropic and homogeneous media:

$\begin{matrix}{{{{\nabla{\times \overset{\rightarrow}{E}}} = {{- j}\; \omega \overset{\rightarrow}{B}}}{{\nabla{\times \overset{\rightarrow}{B}}} = {{\mu \; \overset{\rightarrow}{J}} + {j\; {\omega ɛ\mu}\; \overset{\rightarrow}{E}}}}}{{\nabla{\cdot \overset{\rightarrow}{E}}} = \frac{\rho}{ɛ}}{{\nabla{\cdot B}} = 0}} & (1)\end{matrix}$

where ε, μ, {right arrow over (J)} and ρ are permittivity, permeability,angular frequency, free volumetric electric current and chargedensities, respectively, and j=√{square root over (−1)} (noting that inthe frequency domain of time-harmonic fields,

$\left. \frac{\partial}{\partial t}\rightarrow{j\; \omega} \right).$

This difficulty in direct calculation comes from the measurable EMsources i.e. ρ and {right arrow over (J)} which do not give {right arrowover (E)} and {right arrow over (B)} in a straightforward manner throughthe following decoupled equations for {right arrow over (E)} and {rightarrow over (B)} derived from (1):

$\begin{matrix}{{{{{\nabla^{2}E_{i}} + {k^{2}E_{i}}} = {\left. {{j\; \omega \; \mu \; J_{i}} + \left( {\nabla\left( \frac{\rho}{ɛ} \right)} \right)_{i}}\Rightarrow i \right. = x}},y,z}{{boundary}\mspace{20mu} {conditions}\text{:}\left\{ \begin{matrix}{E_{t}^{-} = E_{t}^{+}} \\{{{ɛ_{r}^{-}E_{n}^{-}} - {ɛ_{r}^{+}E_{n}^{+}}} = \frac{\rho_{s}}{ɛ_{0}}}\end{matrix} \right.}} & \left( {2a} \right) \\{{{{\nabla^{2}B_{i}} + {k^{2}B_{i}}} = {\left. {\frac{- 1}{j\; {\omega ɛ}}\; {\mu \left( {\nabla{\times J}} \right)}_{i}}\Rightarrow i \right. = x}},y,{z{boundary}\mspace{20mu} {conditions}\text{:}\left\{ \begin{matrix}{{\frac{B_{t}^{-}}{\mu_{r}^{-}}\frac{B_{t}^{+}}{\mu_{r}^{+}}} = {\mu_{0}J_{s}}} \\{B_{n}^{-} = B_{n}^{+}}\end{matrix} \right.}} & \left( {2b} \right)\end{matrix}$

where the superscripts “+” and “−” denote the corresponding quantityright after and right before the boundary, respectively, ρ_(s) and J_(s)are the surface charge and current densities, t and n stand for thetangential and normal directions, and k=ω√{square root over (εμ)} is thewave number. The vector derivatives of measurable EM sources that appearin equations (2a) and (2b) make them difficult to solve directly,although they are still solvable at the expense of increasedcomputational time in calculating the vector derivatives. TheFinite-Difference Time-Domain method (“FDTD”), on the other hand,directly deals with Equation (1), utilizing the Yee method to avoid thedifficulties in Equation (2). However, as the FDTD deals with thecoupled equations, it must simultaneously solve 6 field components,three magnetic and three electric ones (instead of the 3 fieldcomponents in Equation (2)).

To overcome the above difficulty, the fields can be indirectlycalculated by defining two ‘auxiliary’ potentials, v and {right arrowover (A)}, as the scalar electric potential and the vector magneticpotential, respectively, forming a set of four variables (A_(x), A_(y),A_(z), v) that can be directly determined via EM measurable quantities(not their vector derivatives) as

$\begin{matrix}{{{{{\nabla^{2}A_{i}} + {k^{2}A_{i}}} = {\left. {\mu \; J_{i}}\Rightarrow i \right. = x}},y,z}{{{\nabla^{2}v} + {k^{2}v}} = {- \frac{\rho}{ɛ}}}} & \left( {3a} \right)\end{matrix}$

and the corresponding boundary conditions as follows (more details aboutthese boundary conditions can be found in Appendix I):

$\begin{matrix}\left\{ \begin{matrix}{{third}\mspace{14mu} {kind}\mspace{14mu} {on}\mspace{14mu} {truncating}\mspace{14mu} {boundaries}\text{:}\left\{ \begin{matrix}{{{n \cdot {\nabla v}} + {qv}} \approx 0} \\{{{n \cdot {\nabla A_{i}}} + {p\; A_{i}}} \approx 0}\end{matrix} \right.} \\{{continuity}\mspace{14mu} {on}\mspace{14mu} {dielectrics}\text{:}\left\{ \begin{matrix}{v^{-} = {{{v^{+}\&}\mspace{20mu} ɛ_{r}^{-}\frac{\partial v^{-}}{\partial n}} = {ɛ_{r}^{+}\frac{\partial v^{+}}{\partial n}}}} \\{{n \times A_{i}^{-}} = {n \times A_{i}^{+}}} \\{{\frac{1}{\mu_{r}^{-}}\left( {\frac{\partial A_{n}^{-}}{\partial i} - \frac{\partial A_{i}^{-}}{\partial n}} \right)} = {\frac{1}{\mu_{r}^{+}}\left( {\frac{\partial A_{n}^{-}}{\partial i} - \frac{\partial A_{i}^{+}}{\partial n}} \right)}}\end{matrix} \right.} \\{{discontinuity}\mspace{14mu} {on}\mspace{14mu} {conductors}\text{:}\left\{ \begin{matrix}{{{ɛ_{r}^{-}\frac{\partial v^{-}}{\partial n}} - {ɛ_{r}^{+}\frac{\partial v^{+}}{\partial n}}} = \frac{\rho_{s}}{ɛ_{0}}} \\{\frac{1}{\mu_{r}^{-}}\left( {\frac{\partial A_{n}^{-}}{\partial i} - \frac{\partial A_{i}^{-}}{\partial n}} \right)} \\{{{- \frac{1}{\mu_{r}^{+}}}\left( {\frac{\partial A_{n}^{-}}{\partial i} - \frac{\partial A_{i}^{+}}{\partial n}} \right)} = {\mu_{0}J_{s}}}\end{matrix} \right.} \\{{Dirichlet}\mspace{14mu} {on}\mspace{14mu} {boundaries}\mspace{14mu} {with}\mspace{14mu} {given}\mspace{14mu} {value}\text{:}\left\{ \begin{matrix}{v = f} \\{A_{i} = g_{i}}\end{matrix} \right.}\end{matrix} \right. & \left( {3b} \right)\end{matrix}$

where ε, and μ, are the relative permittivity and permeability,respectively. The functions f, g, q and p are defined below. Theboundary conditions in (3b) are classified in more detail below (withrespect to the boundary conditions in (2)) for different boundaries. Theelectric field and magnetic flux density are then derived from theauxiliary potentials described above as:

{right arrow over (E)}=−∇v−jω{right arrow over (A)}

{right arrow over (B)}=∇×{right arrow over (A)}  (4)

The boundary conditions in (3b) are fundamental to every EM problem, andare described below.

-   -   The Continuity Boundary Condition: As seen in Eq. (4), the        electric or magnetic potential cannot be discontinuous, because        this would cause their vector derivatives (which are in Eq. (4))        to be infinite, in contradiction with nature. Consequently, the        auxiliary potentials are continuous on boundaries; i.e., v⁻=v⁺,        A_(i) ⁻=A_(i) ⁺.    -   The Dirichlet Boundary Condition: There are some boundaries,        particularly ports and conductors, on which the potential has a        fixed and known value. For example, consider a dipole antenna        disposed along a z-axis and excited by a sinusoidal voltage at        its port. A sinusoidal voltage distribution is thus formed on        the port of the antenna. This boundary condition is known as the        Dirichlet condition, and can be described on the antenna port as        v=f=V₀ sin(kz), where V₀ is a constant. This voltage applies a        surface charge on the antenna

${\frac{\partial v}{\partial n} = \frac{\rho_{s}}{ɛ}},$

while the voltage over the antenna conducting surface (not the port)remains constant; i.e., V₀ sin(kz₀), where z₀ is the point in spacewhere the antenna conductor is connected to the port. This remains validfor {right arrow over (A)} if the antenna is excited by a sinusoidalcurrent instead of voltage, which results in a magnetic potential at theport of A_(z)=g_(z)=A₀ sin(kz), where A₀ is a constant. Thecorresponding surface current on the conducting surface of the antennais

$\frac{\partial A}{\partial n} = {\mu \; {J_{s}.}}$

-   -   (Absorbing) Boundary Conditions of the Third Kind: Absorbing or        truncating boundary conditions are derived from the asymptotic        behavior of potentials at large distances from EM sources. The        absorbing boundary condition is used to truncate the        computational domain and make it computationally finite in terms        of memory demands. This condition can be derived as follows,        beginning with the electric potential. At large distances, the        electric potential approximately propagates as:

$\begin{matrix}{\left. {v \approx {V_{0}\frac{e^{- {jkr}}}{r}}}\rightarrow\left. {\frac{\partial v}{\partial r} \approx {{- {{jk}\left( {V_{0}\frac{e^{- {jkr}}}{r}} \right)}} - {\frac{1}{r}\left( {V_{0}\frac{e^{- {jkr}}}{r}} \right)}}}\rightarrow\left. {\frac{\partial v}{\partial r} \approx {{- {jkv}} - {\frac{1}{r}v}}}\rightarrow{\frac{\partial v}{\partial r} + {\left( {{jk} + \frac{1}{r}} \right)v}} \right. \right. \right. = {\left. {{{n \cdot {\nabla v}} + {qv}} \approx 0}\rightarrow q \right. = {{jk} + \frac{1}{r}}}} & (5)\end{matrix}$

Hence, the absorbing or truncating boundary condition for electricpotential with a specific value for q is derived. The same approach isvalid for magnetic potential (this general type of the boundarycondition is known in the art as the “third kind” boundary condition).Of course, this boundary condition introduces some error, and thereexist more accurate types of absorbing boundary conditions, includingsome with higher order approximations, or perfectly matched layers,which are not discussed herein.

-   -   Boundary Conditions on Normal Derivatives: As a direct        conclusion from the electromagnetic boundary conditions on        {right arrow over (E)} and {right arrow over (B)} in Eq. (2),        the normal derivatives of electric and magnetic potentials obey        a specific condition on boundaries. Since the electromagnetic        boundary conditions on normal derivatives of auxiliary        potentials are rarely mentioned in the literature and the        boundary conditions are mainly described in terms of {right        arrow over (E)} and {right arrow over (B)}, in Appendix I the        proof of the boundary conditions on normal derivatives in Eq.        (3b) are given.

When the partial differential equations of Eq. (3a) are supplemented bythe boundary conditions in Eq. (3b), a unique solution for either v or{right arrow over (A)} is guaranteed.

Having described above the equations of the four-potential definition ofelectromagnetics and their boundary conditions, a computationalelectromagnetics process using only two potentials in accordance withembodiments of the present invention will now be described.

The continuity equation in electromagnetics establishes that, for timeharmonic electromagnetic fields, the current and charge densities arelinked together as follows:

∇·J=−jωρ  (6)

Consequently, if there is no current in the problem domain, the chargewill be zero as well. This ‘source-free’ condition covers an importantclass of electromagnetic problems, as in many microwave engineeringproblems, the surface currents and charge densities are key. Insource-free electromagnetic problems, Eq. (1) can be rewritten as:

∇×{right arrow over (E)}=−jω{right arrow over (B)}

∇×{right arrow over (B)}=jωεμ{right arrow over (E)}

∇·{right arrow over (E)}=0

∇·{right arrow over (B)}=0  (7)

which, in turn, rewrites Eq. (2) as:

$\begin{matrix}{{{{{\nabla^{2}E_{i}} + {k^{2}E_{i}}} = {\left. 0\Rightarrow i \right. = x}},y,z}{{boundary}\mspace{20mu} {conditions}\text{:}\left\{ {{{{\begin{matrix}{E_{t}^{-} = E_{t}^{+}} \\{{{ɛ_{r}^{-}E_{n}^{-}} - {ɛ_{r}^{+}E_{n}^{+}}} = \frac{\rho_{s}}{ɛ_{0}}}\end{matrix}{\nabla^{2}B_{i}}} + {k^{2}B_{i}}} = {\left. 0\Rightarrow i \right. = x}},y,{z{boundary}\mspace{20mu} {conditions}\text{:}\left\{ \begin{matrix}{{\frac{B_{t}^{-}}{\mu_{r}^{-}}\frac{B_{t}^{+}}{\mu_{r}^{+}}} = {\mu_{0}J_{s}}} \\{B_{n}^{-} = B_{n}^{+}}\end{matrix} \right.}} \right.}} & (8)\end{matrix}$

The source-free equation (8) can still be solved through either of thefollowing two methods. In the first method, the fields are indirectlycalculated by solving for the scalar electric and magnetic vectorpotentials A_(x), A_(y), A_(z), v; i.e.:

∇² A _(i) +k ² A _(i)=0⇒i=x, y, z

∇² v+k ² v=0  (9)

whose boundary conditions are given by Eq. (3b). In this case, theequations of the four potentials must be solved, while there is nosource in the problem domain.

Alternatively, the second method involves direct calculation of the mainfields (say {right arrow over (E)}) containing 3 main components in Eq.(8), and then calculating {right arrow over (B)} as

$\overset{\rightarrow}{B} = {\frac{1}{{- j}\; \omega}{\nabla{\times {\overset{\rightarrow}{E}.}}}}$

This second method is the preferred choice, because only the threeequations of the main vector fields need to be solved.

Equation (3a) indicates that {right arrow over (A)} originates fromelectric currents. Although Eq. (9) represents a source-free problem, itnevertheless requires solving the magnetic vector potential for thethree components A_(x), A_(y), A_(z), and consequently the second of theabove methods is preferably used to solve source-free problems becauseit is computationally less demanding.

In view of the requirement to solve the equations of 3 scalar potentials(A_(x), A_(y), A_(z)) in electromagnetic problems originating from{right arrow over (B)}=∇×{right arrow over (A)}, together with the factthat {right arrow over (B)}=∇×{right arrow over (A)} is only one of thepossible definitions for {right arrow over (B)}, the inventorsquestioned whether {right arrow over (B)}=∇×{right arrow over (A)} isalways the best definition for solving EM problems, or whether there maybe cases where it is better (in terms of computation time and possiblyalso the interpretation of the problem) to use other possibledefinitions. The inventors explored alternative definitions for {rightarrow over (B)} and assessed their suitability and practicality forsolving EM problems.

The definition of B in Eq. (4) arises from its divergence-free identity(∇·{right arrow over (B)}=0), which allows {right arrow over (B)} to bedefined as any divergence-free ‘auxiliary’ field, and indeed any vectorfield satisfying the divergence-free property can be considered as asolution for {right arrow over (B)}. Recognizing this fact, and inaccordance with the described embodiments of the present invention, themagnetic field {right arrow over (B)} is defined by two scalarpotentials, as follows:

{right arrow over (B)}=∇φ×∇ψ  (10)

which clearly satisfies the requirement that ∇·{right arrow over (B)}=0,and with straightforward vector calculus can be rewritten as:

$\begin{matrix}{\overset{\rightarrow}{B} = {\frac{1}{2}{\nabla{\times \left( {{\phi {\nabla\psi}} - {\psi {\nabla\phi}}} \right)}}}} & (11)\end{matrix}$

Substituting {right arrow over (B)} from equation (4) into equation (11)and solving for {right arrow over (A)}, yields:

$\begin{matrix}{\overset{\rightarrow}{A} = {\frac{1}{2}\left( {{\phi {\nabla\psi}} - {\psi {\nabla\phi}}} \right)}} & (12)\end{matrix}$

Although this establishes the relationship between (φ, ψ) and {rightarrow over (A)}, in order to rewrite both {right arrow over (E)} and{right arrow over (B)} in Eq. (4) according to this pair of scalarpotentials, their relation with the scalar potential v should bedetermined. To this end, and as the divergence of {right arrow over (A)}can be freely defined, a value for ∇ {right arrow over (A)} thatsimplifies the field calculations is utilized. In classicalelectromagnetics, this value is established according to the followingequation, known as the Lorentz gauge condition:

∇·{right arrow over (A)}=−jωεμv  (13)

The Lorenz gauge condition of Eq. (13) relates the voltage to magneticpotential in an arbitrary manner to simplify the equations. Hence, tofind the relation between the pair scalars and the electric potential,the same procedure can be followed by equating the right-hand side ofEq. (13) to an arbitrary term containing the pair potentials. Thisrelation, while arbitrary, does not violate the electromagneticframework because it still satisfies Maxwell's equations. In thefollowing, two alternative forms of the definition are given to providetwo possible solutions.

Definition #1:

In the first possible definition, the right hand side of (13) is definedas follows:

$\begin{matrix}{{\nabla{\cdot \overset{\rightarrow}{A}}} = {- \frac{\omega^{2}{ɛ\mu\phi\psi}}{2}}} & (14)\end{matrix}$

An important point in Eq. (14) is that there is no change in the Lorenzgauge condition, because Eq. (14) is assumed to be an identicalrepresentation of Eq. (13), but in terms of the pair potentials. Thereason for this specific chosen value is discussed below where itsimplifies the governing equations on pair scalar potentials. It followsfrom equations (13) and (14) that the scalar electric potential v isrelated to the scalar pair-potentials φ and ψ according to:

$\begin{matrix}{v = {{- \frac{j\; \omega}{2}}{\phi\psi}}} & (15)\end{matrix}$

Therefore, the electric field {right arrow over (E)} can be rewritten interms of the pair-potentials by substituting equations (15) and (12)into equation (4) to provide:

{right arrow over (E)}=jωψ∇φ  (16)

Now that equations (10) and (16) express {right arrow over (E)} and{right arrow over (B)} in terms of the pair potentials, the governingequations on the pair potentials can be specified to replace theclassical auxiliary potential equations on {right arrow over (A)} and v.By substituting Eq. (12) in Eq. (14), the governing equation on thescalar potentials is provided as:

$\begin{matrix}{{\frac{1}{\psi}\left( {{\nabla^{2}\psi} + {k^{2}\psi}} \right)} = {\frac{1}{\phi}{\nabla^{2}\phi}}} & (17)\end{matrix}$

where k=ω√{square root over (εμ)} is the well-known wave number.Equation (17) is still a coupled equation for which the individualidentities of φ and ψ are unknown.

In order to solve (17) to find φ and ψ, it must be separated so thateither φ or ψ has their own equations. To decouple these potentials,(17) is equated to γ² as an arbitrary value, providing the followingseparated equations:

∇²ψ+(k ²−γ²)ψ=0

∇²φ−γ²φ=0  (18)

To calculate EM fields, these separate equations are solved individuallyat step 202, and at step 204 their solutions are substituted intoequations (10) and (16) to respectively determine the correspondingmagnetic and electric fields. However, it should be noted that equation(15) represents the electric potential v in terms of the two scalarpotentials φ and ψ. Since the electric potential is a well-knownquantity in terms of the boundary conditions, there is no need tocalculate both ψ and φ directly, and either one of them, say ψ, can besimply calculated after solving the partial differential equations of vand the other one, say φ, determined as

$\begin{matrix}{\phi = {- \frac{2v}{j\; {\omega\psi}}}} & (19)\end{matrix}$

Equation (19) is a simple division, and its computational time is almostnegligible. Therefore, the focus will be only on the first equation in(18). In order to solve the equation of ψ in (18), the value of γ shouldfirst be determined. Assuming isotropic and homogeneous media, thewavenumber of EM waves (k) is a scalar. This fact is utilized inAppendix II to prove that the only possible value for γ is γ=0, whichshows that φ has a static nature (but clearly with non-trivialsolutions, as its Dirichlet boundary conditions are not zero onconductors. Thus, the next step is to find the distribution of ψ.

In order to find the distribution of ψ, its boundary conditions aredefined so that a unique solution is guaranteed. The boundary conditionsfor ψ are derived as follows.

-   -   The Continuity Boundary Condition: As seen in (10), the gradient        of ψ appears in the equation of {right arrow over (B)}. As        {right arrow over (B)} cannot be infinite anywhere, ψ cannot be        discontinuous anywhere. Thus, ψ⁻=ψ⁺ is applies on boundaries.    -   The Third Kind of Boundary Condition on Tangential Derivatives:        from equation (10):

$\begin{matrix}{\overset{\rightarrow}{B} = {{{\nabla\phi} \times {\nabla\psi}} = {{\begin{matrix}\overset{\rightarrow}{t_{1}} & \overset{\rightarrow}{t_{2}} & \overset{\rightarrow}{n} \\\frac{\partial\phi}{\partial t_{1}} & \frac{\partial\phi}{\partial t_{2}} & \frac{\partial\phi}{\partial n} \\\frac{\partial\psi}{\partial t_{1}} & \frac{\partial\psi}{\partial t_{2}} & \frac{\partial\psi}{\partial n}\end{matrix}} = {{\left( {{\frac{\partial\phi}{\partial t_{2}}\frac{\partial\psi}{\partial n}} - {\frac{\partial\phi}{\partial n}\frac{\partial\psi}{\partial t_{2}}}} \right){\overset{\rightarrow}{t}}_{1}} + {\left( {{\frac{\partial\phi}{\partial n}\frac{\partial\psi}{\partial t_{1}}} - {\frac{\partial\phi}{\partial t_{1}}\frac{\partial\psi}{\partial n}}} \right){\overset{\rightarrow}{t}}_{2}} + {\left( {{\frac{\partial\phi}{\partial t_{1}}\frac{\partial\psi}{\partial t_{2}}} - {\frac{\partial\phi}{\partial t_{2}}\frac{\partial\psi}{\partial t_{2}}}} \right)\overset{\rightarrow}{n}}}}}} & (20)\end{matrix}$

-   -   where t₁, t₂ and n stand for the two tangential components and        the normal component in Cartesian coordinates, respectively.        Substituting equation (19) into equation (20) gives:

$\begin{matrix}{\overset{\rightarrow}{B} = {{\frac{2}{j\; {\omega\psi}}\left( {{{- \frac{\partial v}{\partial t_{2}}}\frac{\partial\psi}{\partial n}} + {\frac{\partial v}{\partial n}\frac{\partial\psi}{\partial t_{2}}}} \right){\overset{\rightarrow}{t}}_{1}} + {\frac{2}{j\; {\omega\psi}}\left( {{{- \frac{\partial v}{\partial n}}\frac{\partial\psi}{\partial t_{1}}} + {\frac{\partial v}{\partial t_{1}}\frac{\partial\psi}{\partial n}}} \right){\overset{\rightarrow}{t}}_{2}} + {\frac{2}{j\; {\omega\psi}}\left( {{{- \frac{\partial v}{\partial t_{1}}}\frac{\partial\psi}{\partial t_{2}}} + {\frac{\partial v}{\partial t_{2}}\frac{\partial\psi}{\partial t_{1}}}} \right)\overset{\rightarrow}{n}}}} & (21)\end{matrix}$

-   -   Now, compare this equation with the boundary conditions of the        magnetic flux density in electromagnetics. On conductors, the        tangential component of the magnetic flux density supports a        surface current (e.g., B_(t) ₁ =μJ, in t₁ direction). Also, the        electric voltage is constant across the conductors, resulting in        zero tangential component for its gradient

$\left( {\frac{\partial v}{\partial t_{1}} = {\frac{\partial v}{\partial t_{2}} = 0}} \right).$

These two conditions assign the following tangential boundary conditionfor ψ on the conductors:

$\begin{matrix}{{\frac{\partial\psi}{\partial t_{2}} - {\left( \frac{j\; \omega \; \mu \; J_{s}}{2\frac{\partial v}{\partial n}} \right)\psi}} = 0} & (22)\end{matrix}$

-   -   Moreover, the normal boundary condition of the electric voltage        in equation (3b) means that its normal derivative supports a        surface charge on the conductors, as follows:

$\begin{matrix}{{\frac{\partial\psi}{\partial t_{2}} - {\left( \frac{j\; {\omega ɛ\mu}\; J_{s}}{2\rho_{s}} \right)\psi}} = 0} & (23)\end{matrix}$

-   -   The relationship between the surface current density and charge        density is J_(s)=jωρ_(s). Substituting this relationship into        equation (23) finalizes the third kind of boundary condition of        ψ over the conductors as:

$\begin{matrix}{{\frac{\partial\psi}{\partial t_{2}} + {\frac{k^{2}}{\underset{h}{\underset{}{2}}}\psi}} = 0} & (24)\end{matrix}$

-   -   Are equations (16) and (24) in contradiction on conductors? The        tangential component of the electric field is always zero on        conductors. Therefore, one may incorrectly derive:

{right arrow over (E)} _(t) =jωψ∇ _(t)φ=0  (25)

-   -   on the conductors, which leads us toward either of the following        (incorrect) conditions on conductors:

ψ=0

or

∇_(t)φ=0  (26)

If ψ=0, then equation (19) would cause ϕ to be infinite, incontradiction to the continuity boundary conditions. Conversely, if∇_(t)φ=0, this condition assigns a constant value to both φ=φ₀ and ψ=ψ₀on conductors. In other words, ∇_(t)ψ=0 which violates equation (24).Consequently, equation (25) cannot be correct, and indeed this incorrectapproach is due to the direct use of equation (16). As shown in AppendixI, on (conducting) boundaries, only the electric potential contributesto derive its normal (or tangential) boundary conditions from theelectric fields, and the magnetic potential {right arrow over (A)} doesnot have any contribution (it is canceled in the calculation process).Hence, as equation (16) is based on both {right arrow over (A)} and v,its direct use is wrong.

-   -   The Third Kind or Absorbing Boundary Condition: Since ψ is a        potential function, the same procedure as electric potential is        utilized to approximate ψ at far fields. The behavior of ψ at        far fields is thus approximated by the following third kind of        boundary condition:

$\begin{matrix}{{{n \cdot {\nabla\psi}} + {\underset{\underset{w}{}}{\left( {{jk} + \frac{1}{r}} \right)}\psi}} \approx 0} & (27)\end{matrix}$

-   -   Consequently, a typical absorbing boundary condition for ψ is        derived.    -   Boundary Conditions on Normal Derivatives: Substituting        equation (21) in the tangential boundary conditions of the        magnetic flux density in equation (8), the following relation        between the normal and tangential derivatives of ψ and the        electric voltage on dielectrics is derived:

${\frac{\partial v}{\partial t}\left\lbrack {{\frac{1}{\mu^{+}}\frac{\partial\psi^{+}}{\partial n}} - {\frac{1}{\mu^{-}}\frac{\partial\psi^{-}}{\partial n}}} \right\rbrack} = {{\frac{\partial v}{\partial n}\left\lbrack {{\frac{1}{\mu^{+}}\frac{\partial\psi^{+}}{\partial t}} - {\frac{ɛ^{+}}{ɛ^{-}\mu^{-}}\frac{\partial\psi^{-}}{\partial t}}} \right\rbrack}.}$

-   -   In this equation, the normal and tangential derivatives are        linked together. However, as the normal and tangential        directions are perpendicular and thus independent from each        other, this equation is equated to zero to apply the tangential        and normal derivative boundary conditions, independently.        Experimental observations show that the desired electromagnetic        fields are derived under the above condition. Applying this        condition gives:

${{\frac{1}{\mu^{+}}\frac{\partial\psi^{+}}{\partial n}} = {\frac{1}{\mu^{-}}\frac{\partial\psi^{-}}{\partial n}}},{{\frac{1}{\mu^{+}ɛ^{+}}\frac{\partial\psi^{+}}{\partial t}} = {\frac{1}{ɛ^{-}\mu^{-}}{\frac{\partial\psi^{-}}{\partial t}.}}}$

-   -   The above boundary condition thus describes the normal and        tangential derivative components of ψ over the dielectric        interfaces.

Having derived the boundary conditions for ψ as a new potential, itscorresponding equation will provide a unique solution. Consequently, bydetermining the pair scalar potentials through solving the governingequations on v and ψ supplemented by their boundary conditions, themagnetic flux density and electric field in source-free problems can becalculated according to equations (10) and (16), respectively, usingonly two scalar potentials (instead of the 4 auxiliary or 3 maincomponents of the prior art).

By way of summary, the governing equations on v and ψ are:

∇² ψ+k ²ψ=0

∇² v+k ² v=0  (28a)

Since the boundary conditions of the electric potential have alreadybeen summarized above in equation (3b), only the boundary conditions ofψ are summarized here as:

$\begin{matrix}\left\{ \begin{matrix}{{{{third}\mspace{14mu} {kind}\mspace{14mu} {on}\mspace{14mu} {truncating}\mspace{14mu} {boundaries}\text{:}\mspace{14mu} {n \cdot {\nabla\psi}}} + {w\; \psi}} \approx 0} \\{{continuity}\mspace{14mu} {on}\mspace{14mu} {dielectrics}\text{:}\mspace{14mu} \left\{ \begin{matrix}{{\frac{1}{\mu_{r}^{-}}\frac{\partial\psi^{-}}{\partial n}} = {\frac{1}{\mu_{r}^{+}}\frac{\partial\psi^{+}}{\partial n}}} \\{\psi^{-} = \psi^{+}} \\{{\frac{1}{\mu_{r}^{+}ɛ_{r}^{+}}\frac{\partial\psi^{+}}{\partial t}} = {\frac{1}{ɛ_{r}^{-}\mu_{r}^{-}}\frac{\partial\psi^{-}}{\partial t}}}\end{matrix} \right.} \\{{{{third}\mspace{14mu} {kind}\mspace{14mu} {on}\mspace{11mu} {conductors}\text{:}\mspace{14mu} n \times {\nabla\psi}} + {h\; \psi}} = 0}\end{matrix} \right. & \left( {28b} \right)\end{matrix}$

Definition #2:

The second possible definition is given as:

∇·{right arrow over (A)}=−ω ²εμφψ  (29)

An important point in equation (29) is that there is no change in theLorenz gauge condition because equation (29) is assumed to be anidentical representation of equation (13), but in terms of the pairpotentials. This specific form is chosen to simplify the governingequations on the pair scalar potentials, as described below.

Equations (13) and (29) can be used to find the relation between thescalar electric potential v and the scalar pair-potentials φ and ψ as:

v=−jωφψ  (30)

Therefore, the electric field {right arrow over (E)} can now berewritten according to the pair-potentials by substituting equations(30) and (12) into equation (4) as:

$\begin{matrix}{\overset{\rightharpoonup}{E} = {\frac{j\; \omega}{2}\left( {{\phi {\nabla\psi}} + {3\psi {\nabla\phi}}} \right)}} & (31)\end{matrix}$

Thus equations (10) and (31) respectively define the magnetic andelectric fields in terms of the pair-scalar potentials. With thesedefinitions, the governing equations on the presented potentials can bespecified so that they can replace the classical auxiliary potentialequations on {right arrow over (A)} and v. By substituting equation (12)in equation (29), the governing equation on the presented potentials isderived as:

$\begin{matrix}{{\frac{1}{\psi}\left( {{\nabla^{2}\psi} + {2k^{2}\psi}} \right)} = {\frac{1}{\phi}{\nabla^{2}\phi}}} & (32)\end{matrix}$

However, equation (32) is still a coupled equation for φ and ψ. Todecouple these potentials, equation (32) is equated to γ² as anarbitrary function leading us towards the following separated equations:

∇²ψ+(2k ²−γ²)ψ=0

∇²φ−γ²φ=0  (33)

In practice, these separate equations are solved individually and theirsolutions substituted into equations (10) and (31) to find the magneticand electric fields, respectively. However, equation (30) defines theelectric potential v as another scalar in terms of φ and ψ. Since theelectric potential v is a well-known quantity in terms of sources andboundary conditions, there is no need to calculate both ψ and φdirectly. One of them, say φ, can be simply calculated after solving thepartial differential equations of v and the other one, say ψ, determinedas:

$\begin{matrix}{\phi = {- \frac{v}{j\; {\omega\psi}}}} & (34)\end{matrix}$

Clearly, equation (34) is a simple division, and its computational timeis therefore almost negligible. Therefore, the focus will be only on thefirst equation in (33).

The first equation in (33) is an eigenvalue equation with wavenumber√{square root over (2k)} and eigenvalue γ. The meaning of this equationis discussed later, after describing the boundary conditions of ψ. Thecorresponding distribution of ψ for the eigenvalue is called theeigenstate of ψ. In every eigenvalue problem, the first step is to findthe eigenvalue itself. Then, the corresponding eigenstates of thiseigenvalue are found. Assuming isotropic and homogeneous media, thewavenumber of EM waves (k) is therefore a scalar. This fact is utilizedin the Appendix IV to prove that the only possible value for theeigenvalue is

$\gamma = {{- \frac{j}{2}}{k.}}$

Interestingly, if the eigenvalues of ψ are independently derived usingFEM, not using the simple approach given in the Appendix IV, thedominant eigenvalue will be

${\gamma = {{- \frac{j}{2}}k}},$

indicating that the two approaches coincide and agree that the higherorder eigenvalues are not applicable. Thus, it is only necessary to findthe eigenstates of ψ, because the eigenvalue is already known.

To find the eigenstates of ψ, its boundary conditions are defined toguarantee a unique solution. The boundary conditions for it are derivedas follows:

-   -   The Continuity Condition: As seen in equations (10) and (31),        the gradient of it appears in the equations for {right arrow        over (E)} and {right arrow over (B)}. As {right arrow over (E)}        and {right arrow over (B)} cannot be infinite anywhere, ψ cannot        be discontinuous anywhere. Thus, ψ⁻=ψ⁺ on boundaries.    -   The Dirichlet Boundary Condition: In eigenvalue problems, there        is no excitation in the problem domain. In other words, the        exciting ports are all off. This assigns a zero value ψ it on        conductors. The Dirichlet boundary condition on conductors for ψ        is thus ψ=h=0 (where h is a symbolic function).    -   The Absorbing Boundary Condition: The same procedure as electric        potential is utilized to approximate ψ at far fields. Obviously,        the eigenstates of ψ are standing waves as they are derived from        an eigenvalue problem without any excitation term. Since the        standing waves do not decay at far fields, their behaviour at        far fields is approximated by:

$\left. {\psi \approx {\psi_{0}e^{{- j}\sqrt{{2k^{2}} - \gamma^{2}}r}}}\rightarrow\left. {\frac{\partial\psi}{\partial n} \approx {{- j}\sqrt{{2k^{2}} - \gamma^{2}}\psi}}\rightarrow{{{n \cdot {\nabla\psi}} + {\underset{\underset{w}{}}{j\sqrt{{2k^{2}} - \gamma^{2}}}\psi}} \approx 0} \right. \right.$

-   -   Consequently, a typical absorbing boundary condition for ψ is        derived.    -   Boundary Conditions on Normal and Tangential Derivatives:        Substituting (12) in the normal boundary conditions of the        magnetic vector potential in (3b), and then, substituting for φ        using (34), the following relation between the normal and        tangential derivatives of ψ and voltage on dielectrics is        derived:

$\left\lbrack {{\frac{1}{\mu^{+}}\frac{\partial\psi^{+}}{\partial n}} - {\frac{1}{\mu^{-}}\frac{\partial\psi^{-}}{\partial n}}} \right\rbrack = {{\frac{\frac{\partial v}{\partial n}}{\frac{\partial v}{\partial t}}\left\lbrack {{\frac{1}{\mu^{+}}\frac{\partial\psi^{+}}{\partial t}} - {\frac{ɛ^{+}}{ɛ^{-}\mu^{-}}\frac{\partial\psi^{-}}{\partial t}}} \right\rbrack}.}$

-   -   In this equation, the normal and tangential derivatives are        linked together. However, as the normal and tangential        directions are perpendicular and thus independent from each        other, the above equation is equated to zero to apply the        tangential and normal derivative boundary conditions,        separately. Experimental observations show that the desired        electromagnetic fields are derived under the above separation.        Applying this condition gives:

${{\frac{1}{\mu^{+}}\frac{\partial\psi^{+}}{\partial n}} = {\frac{1}{\mu^{-}}\frac{\partial\psi^{-}}{\partial n}}},{{\frac{1}{\mu^{+}ɛ^{+}}\frac{\partial\psi^{+}}{\partial t}} = {\frac{1}{ɛ^{-}\mu^{-}}{\frac{\partial\psi^{-}}{\partial t}.}}}$

-   -   The above boundary condition thus describes the normal and        tangential derivative components of ψ over the dielectric        interfaces.

The boundary conditions have now been derived for ψ, and hence itscorresponding equation will give a unique solution. Because the equationand boundary conditions of ψ are reminiscent of the properties andbehaviour of energy density, this is the physical quantity allocated toψ in this specification.

By way of summary, in order to solve electromagnetic problems, the firststep is to determine the pair scalar potentials by solving the governingequations on v and ψ supplemented by their boundary conditions. Themagnetic flux density and electric fields are then calculated accordingto equations (10) and (31), respectively, using only two scalars. Thegoverning equations and boundary conditions on v and ψ are:

$\begin{matrix}{{{{\nabla^{2}\psi} + {\left( {{2k^{2}} - \gamma^{2}} \right)\psi}} = 0}{{{\nabla^{2}v} + {k^{2}v}} = {- \frac{\rho + \rho_{J}}{ɛ}}}} & \left( {35a} \right) \\\left\{ \begin{matrix}{{continuity}\mspace{14mu} {on}\mspace{14mu} {dielectrics}\text{:}\mspace{14mu} \left\{ \begin{matrix}{v^{-} = {{{v^{+}\&}\mspace{14mu} ɛ_{r}^{-}\frac{\partial v^{-}}{\partial n}} = {ɛ_{r}^{+}\frac{\partial v^{+}}{\partial n}}}} \\{{\psi^{-} = \psi^{+}}\mspace{211mu}} \\{{{\frac{1}{\mu_{r}^{-}}\frac{\partial\psi^{-}}{\partial n}} = {\frac{1}{\mu_{r}^{+}}\frac{\partial\psi^{+}}{\partial n}}}} \\{{{\frac{1}{\mu_{r}^{+}ɛ_{r}^{+}}\frac{\partial\psi^{+}}{\partial t}} = {\frac{1}{ɛ_{r}^{-}\mu_{r}^{-}}\frac{\partial\psi^{-}}{\partial t}}}\mspace{34mu}}\end{matrix} \right.} \\{{continuity}\mspace{14mu} {on}\mspace{14mu} {conductors}\text{:}\mspace{14mu} \left\{ \begin{matrix}{{v^{-} = v^{+}}\mspace{11mu}} \\{\psi^{-} = \psi^{+}}\end{matrix} \right.} \\{{Dirichlet}\mspace{14mu} {on}\mspace{14mu} {boundaries}\mspace{14mu} {with}\mspace{14mu} {given}\mspace{14mu} {value}\text{:}\mspace{14mu} \left\{ \begin{matrix}{v = f} \\{\psi = h}\end{matrix} \right.} \\{{absorbing}\mspace{14mu} {on}\mspace{14mu} {truncating}\mspace{14mu} {boundaries}\text{:}\mspace{14mu} \left\{ \begin{matrix}{{{{n \cdot {\nabla v}} + {qv}} \approx 0}\mspace{11mu}} \\{{{n \cdot {\nabla\psi}} + {w\; \psi}} \approx 0}\end{matrix} \right.}\end{matrix} \right. & \left( {35b} \right)\end{matrix}$

The change in the right-hand side of the equation of voltage i.e. ρ_(J),is discussed further. As seen in the derived equations (35a), while thecharge density contributes in the 2-scalar equations as well as theclassical ones, the current density J has disappeared. However, thisdisappearance does not mean ignoring the effect of current densities.The continuity equation in electromagnetics, i.e.:

∇·J=−jωρ _(J)  (36)

tells us that the current and charge densities for time harmonic fieldsare linked together, and knowing one of them means the other can bedetermined. Therefore, if there is any current in the problem domain,this current is identical to an equivalent charge of

$\rho_{J} = {{- \frac{1}{j\; \omega}}{\nabla{\cdot J}}}$

and this charge and other free charges are considered in the second ofequations (35a).

In practice, different computational methods can be used to solve theequations above, but the most common of these methods are finitedifference and finite element methods. Although these methods arestandard methods well known to those skilled in the art, and may beprovided, for example, by existing software packages (as they are in theExamples described below), the application of each of these methods tosolve electromagnetic problems using the equations described above isbriefly described below to illustrate how they are applied in practice.Appendix III discusses the FEM implementation and formulation fornumerically solving the electric potential equation and the equation ofψ, respectively.

In some embodiments, as shown in FIG. 1, a computationalelectromagnetics system 100 includes an Intel Architecture computer ordata processing system, and executes a computational electromagneticsprocess as shown in FIG. 2 and described herein, the process beingimplemented in the form of programming instructions of one or moresoftware modules 102 stored on non-volatile (e.g., hard disk orsolid-state drive) storage 104 of or associated with the computer system100, as shown in FIG. 1, and configured to solve the equations describedabove, namely either:

-   -   (i) equations (18) or (28a) subject to boundary conditions of        (3b) and (28b) (if using definition #1); or    -   (ii) equations (35a) subject to boundary conditions of (3b) and        (35b) (if using definition #2),

and in either case then solving equations (10) and (16). However, itwill be apparent that the computational electromagnetics process couldalternatively be implemented in the form of one or more dedicatedhardware components, such as application-specific integrated circuits(ASICs) and/or as configuration data of at least one field programmablegate array (FPGA), for example.

In the embodiment of FIG. 1, the system 100 includes a random accessmemory (RAM) 106, at least one processor 108, and external interfaces110, 112, 114, all interconnected by a bus 116. The external interfacesmay include universal serial bus (USB) interfaces 110, at least one ofwhich is connected to a keyboard 118 and a pointing device such as amouse 119, a network interface connector (NIC) 112 which connects thesystem 100 to a communications network such as the Internet 120, and adisplay adapter 114, which is connected to a display device such as anLCD panel display 122. The system 100 may also include one or morestandard software modules 1 to 130, including an operating system 124such as Linux or Microsoft Windows, and optionally web server software126 such as Apache, available at http://www.apache.org. In someembodiments, the computational electromagnetics process described hereincan be embedded in a software package for solving generalelectromagnetics problems, such as any of the state of the artcommercial software packages known to those skilled in the art,including HFSS, CST MS, ADS, COMSOL, XFdtd and OmniSim FDTD, which asdescribed above implement one or more of FDM, FEM, MOM, and meshlessmethods known to those skilled in the art for solving the equationsdescribed herein.

The computational electromagnetics process and system described hereinenable the simulation of electric and/or magnetic fields in a fractionof the time taken by prior art systems and processes to generate resultsto the same accuracy and using the same computational hardware. Forexample, one particular simulation of the safety of an airplane withrespect to electromagnetic fields currently takes several days tocalculate using existing state of the art methods, whereas the samesimulation can be performed in less than a day by the computationalelectromagnetics system and process described herein. In the biomedicalapplication described below, the simulation described below took severalhours using prior art methods, but less than an hour using thecomputational electromagnetics system and process described herein.Moreover, the improved efficiency of calculation can be even moresignificant in optimization simulations that may require, for example,hundreds of iterations to complete. For example, a hyperthermia cancertreatment simulation requiring hundreds of iterations of electromagneticcalculations on a patient-specific breast model (to guarantee properbeamforming and focus the energy at the tumor without any hot spots inhealthy tissue) took more than a week to calculate using prior artmethods on a general purpose personal computer, but was completed inless than a day by the computational electromagnetics process describedherein and using the same personal computer. The computationalelectromagnetics process and system described herein is generallyapplicable to any source-free problem and to problems with sourceshaving known surface current distributions (which are responsible forgenerating the electromagnetic fields). Methods for predicting surfacecurrent distributions in electromagnetic problems with sources are knownto those skilled in the art, including the method described above in thecontext of boundary conditions.

EXAMPLES

In order to validate and demonstrate the power of the described methods,several electromagnetic problems were simulated by prior art classicalmethods and by the new computational electromagnetics process describedherein, using the commercial Finite Element Modelling (FEM) softwarepackage COMSOL Multiphysics®, as described athttps://www.consol.com/comsol-multiphysics.

The COMSOL software package was chosen because it not only providescomponents configured to solve pre-defined problems in electromagneticsand other subject areas using state of the art classical methods, butalso provides an equation-based Mathematics Module modelling componentfor custom simulations, allowing the powerful COMSOL solvers to beapplied to solve essentially any user-defined differential equations andwith any kind of boundary conditions in either 1, 2, or 3D domains. (Incontrast, other commercial software packages, such as ADS, CST and HFSSare closed packages that do not allow users to solve arbitraryuser-defined problems.)

In the following examples, each electromagnetic problem was solved firstby classical methods using the Electromagnetic Waves component of COMSOLMultiphysics®, and then by the method described herein using the“Equation-based Modelling” component of COMSOL Multiphysics®. Becausequalitative comparisons between these two formulations are difficult(e.g., by comparing the color bars in the figures), the followingdifference definition was used to quantitatively demonstrate theaccuracy of the two-scalar method described herein with respect to theclassical solution:

$\begin{matrix}{{Error} = {\frac{\sum\limits_{m}{{{\overset{\rightarrow}{E}}_{classical}^{m} - {\overset{\rightarrow}{E}}_{2 - {variable}}^{m}}}}{M{\sum\limits_{m}{{\overset{\rightarrow}{E}}_{classical}^{m}}}} + \frac{\sum\limits_{m}{{{\overset{\rightarrow}{B}}_{classical}^{m} - {\overset{\rightarrow}{B}}_{2 - {variable}}^{m}}}}{M{\sum\limits_{m}{{\overset{\rightarrow}{B}}_{classical}^{m}}}}}} & (29)\end{matrix}$

where m is the element number of FEM discretization, and M is the totalnumber of elements. It should be borne in mind that the classicalimplementation in COMSOL benefits from physics-based meshing andassociated algorithms that have been optimized to tradeoff betweenaccuracy and computational time, whereas the Mathematics Module used forthe two-scalar method described herein does not benefit from theseoptimizations. Consequently, the reductions in computational timedescribed below are effectively under estimates of the improvements thatare achievable in practice.

To derive electromagnetic fields using the two-scalar method, theequations in (28) and (35), which give identical results, were solved,whereas COMSOL uses the equations in (2) for the classical formulation.All of the simulations were generated by a computer having an IntelCore™ i7 processor operating at 3.60 GHz, and with 16 GB of RAM.

Example I

Dipole Antenna

As shown in FIG. 3(a), a dipole antenna with arm dimension 7 cm×0.5 cmand a port size 0.5 cm×0.5 cm was operated in free space at a frequencyof 1 GHz. The port impedance was 75 Ω in all the simulations, and theantenna surface was considered to be a perfect conductor.

FIGS. 3(b) and 3(c) respectively show the two simulated potentials v andψ (which define the main fields) in the xy, yz, and zx planes. FIGS.3(d) and 3(e) show (in the same three planes) the corresponding electricfields calculated by the classical approach (FIG. 3(d)) and by the2-scalar formulation (FIG. 3(e)). The simulated magnetic flux densitycalculated by the classical method and by the two-scalar method areshown in FIGS. 3(f) and 3(g), respectively.

It is immediately apparent that the fields calculated by thepair-potential method agree very well with those obtained using theclassical method (a quantitative comparison is described below). Slightdifferences are expected due to some physics-based optimizations of theRF Module that do not exist in the Mathematics Module.

Example II

Dipole Antenna With Scatterer

In a second example, a conducting sphere with a radius r=2 cm was placedin front of the dipole antenna of Example I, as shown in FIG. 4(a). Theresulting scattering effects were simulated by both the classical andtwo-scalar methods. FIGS. 4(b) and 4(c) respectively show the twosimulated potentials v and yr in the xy, yz, and zx planes. FIGS. 4(d)and 4€ show (in the same three planes) the corresponding electric fieldscalculated by the classical approach (FIG. 4(d)) and by the two-scalarformulation (FIG. 4(e)). The simulated magnetic flux density calculatedby the classical method and by the two-scalar method are shown in FIGS.4(f) and 4(g), respectively.

As with the first Example, the field distributions calculated by the twodifferent methods are in good agreement, and once again the minordifferences may originate from the different optimizations in the RF andMathematics modules in COMSOL, as described above.

Example III

Array of Dipole Antennas Around a Head Model

In a third Example, the dipole antenna described above was replicated toconstruct an elliptic array of 8 identical antennae. As shown in FIG.5(a), this array was arranged to surround a realistically size andshaped model of a human head with average dielectric propertiesε_(r)=44.45, σ=0.94 (S/m). The head model includes the scatterer fromthe previous example as a conducting object embedded within the head andnear the back side of the head. The purpose of this example was todemonstrate the ability of the two-scalar potential method to solvingmultilayer and complex geometry problems. All of the antennassimultaneously illuminated the head at a frequency of 1 GHz. FIGS. 5(b)to 5(g) show the spatial distributions of the scalar potentials, andFIGS. 5(d) to 5(g) show the spatial distributions of the electric fields(FIGS. 5(d) and (e)) and the magnetic flux density (FIGS. 5(f) and (g)),as calculated by the classical and two-scalar potential methods. Atleast qualitatively, the calculated spatial distributions of theelectric and magnetic fields are clearly in agreement between the twomethods.

A. Comparison

Equation (29) was used to quantify the minor differences between theresults calculated by the classical and two-scalar methods, and theresults are summarized in Table I below. It is apparent that thediscrepancies are quite small, and it is considered that these can beattributed to the different optimization algorithms used in thedifferent modules of COMSOL used by the two methods. For example, theinventors have experienced similarly minor discrepancies in solving agiven problem using different computational software packages; e.g.,HFSS, CST, COMSOL and ADS. The differences listed in Table I can bereduced to below 0.3×10⁻³, 2.61×10⁻³, and 3.21×10⁻³ for the dipole,dipole scatterer and head imaging examples, respectively, when usingproper values for the truncating boundary radius, element discretizationand number of meshes (for example increasing radius of the truncatingsphere from 25 cm to 45 cm, using cubic instead of quadratic elementdiscretization or using fine instead of normal meshing).

B. Computational Time Comparison

The computational times for the above examples using both the classicaland 2-scalar methods are summarized in Table II below. As expected,reducing the number of potentials directly impacts the computationaltime. In every case the two-scalar potential method requires lesscomputational time because it deals with solving two potential equationsin the domain instead of solving the three equations of the main vectorfields. This reduction in computational time is more significant inlarge-scale problems where the size of the matrix representation of thesystem of equations becomes larger (See Appendix III for more details).If the matrix size is M×M, the inversion time will be of the third orderi.e. M³ (while some optimizations may reduce this order to M^(2.37)).

TABLE I QUANTITATIVE COMPARISON FOR EXAMPLES IN SECTION IV USING THECLASSICAL AND 2-SCALAR APPROACHES dipole and dipole antenna scattererhead imaging Difference 0.4 × 10⁻³ 8.1 × 10⁻³ 1.73 × 10⁻² Thesimulations in COMSOL Multiphysics include 3 × 10⁵ linear tetrahedralelements. The error definition is based on (29).

TABLE II COMPUTATIONAL TIME COMPARISON FOR EXAMPLES IN SECTION IV USINGTHE CLASSICAL AND 2-SCALAR APPROACHES dipole and formulation dipoleantenna scatterer head imaging classical 20 sec 31 sec 182 sec 2-scalar18 sec 26 sec 108 sec The simulations in COMSOL Multiphysics include 3 ×10⁵ linear tetrahedral elements.

Many modifications will be apparent to those skilled in the art withoutdeparting from the scope of the present invention.

APPENDICES Appendix I: Normal Boundary Conditions

In order to prove the governing normal boundary conditions on scalarelectric potential seen in (3b), two ways exist. The first one is togeometrically prove them applying the theorems in the calculus ofdifferential and integral, and the second approach is to directlysubstitute the auxiliary potentials in the boundary conditions on themagnetic fields. To show the both methods, we derive the normal boundaryconditions of v through the geometrical way, and that of A through thedirect substitution.

To prove the normal derivative boundary condition of v, a small pillboxwith its top face in medium 1 and its bottom face in medium 2 isconsidered as seen in FIG. 6. The upper and lower faces have an area ΔS,the height of this pillbox is Δh and its volume is V. Taking thevolumetric integrals from the third Maxwell equation on the pillboxyields

$\begin{matrix}{{\nabla{\cdot D}} = {\left. \rho \rightarrow{\int\limits_{V}{\left( {\nabla{\cdot D}} \right){dV}}} \right. = {\left. {\int\limits_{V}{(\rho){dV}}}\rightarrow{ɛ{\int\limits_{V}{\left( {\nabla{\cdot E}} \right){dV}}}} \right. = {\left. {\int\limits_{V}{(\rho){dV}}}\rightarrow{ɛ{\int\limits_{V}{\left( {\nabla{\cdot \left( {{- {\nabla v}} - {j\; \omega \; A}} \right)}} \right){dV}}}} \right. = {\left. {\int\limits_{V}{(\rho){dV}}}\rightarrow{{{- ɛ}{\int\limits_{V}{{\nabla{\cdot \left( {\nabla v} \right)}}{dV}}}} - {j\; {\omega ɛ}{\int\limits_{V}{\nabla{\cdot {AdV}}}}}} \right. = {\int\limits_{V}{(\rho){dV}}}}}}}} & \left( {I{.1}} \right)\end{matrix}$

Using the Lorenz gauge condition in (13), we have

$\begin{matrix}{{{{- ɛ}{\int\limits_{V}{{\nabla{\cdot \left( {\nabla v} \right)}}{dV}}}} - {\omega^{2}ɛ^{2}\mu {\int\limits_{V}{vdV}}}} = {\int\limits_{V}{(\rho){dV}}}} & \left( {I{.2}} \right)\end{matrix}$

To study the behavior of the voltage near the boundaries, we let Δh tobe vanishingly small. Thus, the second integral in the left-hand sidewill be a volume integral on a continuous function (v) for which thevolume is going to zero (v→0). This integral is therefore zero. Theremained terms of (I.2) are

$\begin{matrix}{{{- ɛ}{\int\limits_{V\rightarrow 0}{{\nabla{\cdot \left( {\nabla v} \right)}}{dV}}}} = {\int\limits_{V\rightarrow 0}{(\rho){dV}}}} & \left( {I{.3}} \right)\end{matrix}$

Applying the divergence theorem, the volumetric integral in theleft-hand side is divided into three surface integrals

$\begin{matrix}{{{{- ɛ}\int\limits_{S_{1}}{{\nabla v}.d}{\overset{\rightarrow}{S}}_{1}} - {ɛ{\int\limits_{S_{2}}{{{\nabla v}.d}{\overset{\rightarrow}{S}}_{2}}}} - {ɛ{\int\limits_{S_{3}}{{{\nabla v}.d}{\overset{\rightarrow}{S}}_{3}}}}} = {\int\limits_{V\rightarrow 0}{(\rho){dV}}}} & \left( {I{.4}} \right)\end{matrix}$

where dS₁, dS₂, and dS₃ are the differential of the upper face, lowerface and side face of the pillbox, respectively. Clearly, when Δh→0, theside integral vanishes. Also, it is obvious that dS₁=−dS₂=dS. Hence,

$\begin{matrix}{{{ɛ_{1}{\int\limits_{S}{\left( {{\nabla v}.{\overset{\rightarrow}{a}}_{n}} \right){dS}}}} - {ɛ_{2}{\int\limits_{S}{\left( {{\nabla v}.{\overset{\rightarrow}{a}}_{n}} \right){dS}}}}} = {\int\limits_{V\rightarrow 0}{(\rho){dV}}}} & \left( {I{.5}} \right)\end{matrix}$

where {right arrow over (a)}_(n) is the normal vectors of dS₁. As thepillbox is small, (I.5) is identical to

$\begin{matrix}{{{ɛ_{1}\frac{\partial v_{1}}{\partial n}\Delta \; S} - {ɛ_{2}\frac{\partial v_{2}}{\partial n}\Delta \; S}} = {\left. {\int\limits_{V\rightarrow 0}{(\rho){dV}}}\rightarrow{{ɛ_{1}\frac{\partial v_{1}}{\partial n}} - {ɛ_{2}\frac{\partial v_{2}}{\partial n}}} \right. = {\frac{\int\limits_{V\rightarrow 0}{(\rho){dV}}}{\Delta \; S} = {\left. \rho_{s}\rightarrow{{ɛ_{r}^{-}\frac{\partial v^{-}}{\partial n}} - {ɛ_{r}^{+}\frac{\partial v^{+}}{\partial n}}} \right. = \frac{\rho_{s}}{ɛ_{0}}}}}} & \left( {I{.6}} \right)\end{matrix}$

which finishes the proof. Since ρ_(s)=0 on dielectrics, (I.6) reduces tothe continuity condition in (3b).

In order to prove the boundary condition on magnetic vector potential{right arrow over (A)} in equation (3b), it is more straightforward todirectly deal with the governing boundary condition on the tangentialcomponent of the magnetic field i.e.

{right arrow over (a)} _(n)×(H ₁ −H ₂)=J _(s)  (I.7)

Substituting the definition of B from equation (4) into equation (I.7),we get

$\begin{matrix}{{{\overset{\rightarrow}{a}}_{n} \times \left( {\frac{\nabla{\times A_{1}}}{\mu_{1}} - \frac{\nabla{\times A_{2}}}{\mu_{2}}} \right)} = J_{s}} & \left( {I{.8}} \right)\end{matrix}$

Applying the BAC-CAB rule in vector calculus, equation (I.8) becomes

$\begin{matrix}{\left( {\frac{{\nabla\left( {A_{1}.{\overset{\rightarrow}{a}}_{n}} \right)} - {\nabla{.{{\overset{\rightarrow}{a}}_{n}\left( A_{1} \right)}}}}{\mu_{1}} - \frac{{\nabla\left( {A_{2}.{\overset{\rightarrow}{a}}_{n}} \right)} - {\nabla{.{{\overset{\rightarrow}{a}}_{2}\left( A_{2} \right)}}}}{\mu_{2}}} \right) = J_{s}} & \left( {I{.9}} \right)\end{matrix}$

Therefore,

$\begin{matrix}{{{\frac{1}{\mu_{r}^{-}}\left( {\frac{\partial A_{n}^{-}}{\partial i} - \frac{\partial A_{i}^{-}}{\partial n}} \right)} - {\frac{1}{\mu_{r}^{+}}\left( {\frac{\partial A_{n}^{-}}{\partial i} - \frac{\partial A_{i}^{+}}{\partial n}} \right)}} = {\mu_{0}J_{s}}} & \left( {I{.10}} \right)\end{matrix}$

which finishes the proof. On the dielectric surfaces where J_(s)=0,(I.10) reduces to the continuity condition in equation (3b).

Appendix II: Determining the Value of γ

In order to find γ, the general solution of Equation (18) in Cartesiancoordinates is considered (the same approach is valid for the otherorthogonal coordinate systems). According to the standard method knownto those skilled in the art as “separation of variables”, either φ or ψis written as the product of three functions of one coordinate each (thecoordinate variables are “separated”) as:

ψ(x, y, z)=Ψ_(x)(x)Ψ_(y)(y)Ψ_(z)(z)

φ(x, y, z)=Φ_(x)(x)Φ_(y)(y)Φ_(z)(z)  (II.1)

Substituting (I.1) in Equations (18) yields

$\begin{matrix}\left\{ {\begin{matrix}{{\frac{d^{2}\Psi_{x}}{{dx}^{2}} + {\left( {k_{x}^{2} + \gamma_{x}^{2}} \right)\Psi_{x}}} = 0} \\{{\frac{d^{2}\Psi_{x}}{{dy}^{2}} + {\left( {k_{y}^{2} + \gamma_{y}^{2}} \right)\Psi_{y}}} = 0} \\{{\frac{d^{2}\Psi_{z}}{{dz}^{2}} + {\left( {k_{z}^{2} + \gamma_{z}^{2}} \right)\Psi_{z}}} = 0}\end{matrix},\left\{ \begin{matrix}{{\frac{d^{2}\Phi_{x}}{{dx}^{2}} + {\gamma_{x}^{2}\Phi_{x}}} = 0} \\{{\frac{d^{2}\Phi_{y}}{{dy}^{2}} + {\gamma_{y}^{2}\Phi_{y}}} = 0} \\{{\frac{d^{2}\Phi_{z}}{{dz}^{2}} + {\gamma_{z}^{2}\Phi_{z}}} = 0}\end{matrix} \right.} \right. & \left( {{II}{.2}} \right)\end{matrix}$

where the following relation between the separation parameters ((k_(x),γ_(x)), (k_(y), γ_(y)), (k_(z), γ_(z))) must hold:

$\begin{matrix}{{\left( {k_{x}^{2} + \gamma_{x}^{2}} \right) + \left( {k_{y}^{2} + \gamma_{y}^{2}} \right) + \left( {k_{z}^{2} + \gamma_{z}^{2}} \right)} = \left. {k^{2} + \gamma^{2}}\rightarrow\left\{ \begin{matrix}{{k_{x}^{2} + k_{y}^{2} + k_{z}^{2}} = k^{2}} \\{{\gamma_{x}^{2} + \gamma_{y}^{2} + \gamma_{z}^{2}} = \gamma^{2}}\end{matrix} \right. \right.} & \left( {{II}{.3}} \right)\end{matrix}$

As (II.2) shows, the separated equations are in the same form. Thus,they are usually called “harmonic equations”. In general, the solutionto each of these harmonic equations is expressed by harmonic functionsh_(Ψ), h_(Φ), as:

$\begin{matrix}\left\{ {\begin{matrix}{\Psi_{x} = {h_{\Psi}\left( {\sqrt{k_{x}^{2} + \gamma_{x}^{2}}x} \right)}} \\{\Psi_{y} = {h_{\Psi}\left( {\sqrt{k_{y}^{2} + \gamma_{y}^{2}}y} \right)}} \\{\Psi_{z} = {h_{\Psi}\left( {\sqrt{k_{z}^{2} + \gamma_{z}^{2}}z} \right)}}\end{matrix},\left\{ \begin{matrix}{\Phi_{x} = {h_{\Phi}\left( {\gamma_{x}x} \right)}} \\{\Phi_{y} = {h_{\Phi}\left( {\gamma_{y}y} \right)}} \\{\Phi_{z} = {h_{\Phi}\left( {\gamma_{z}z} \right)}}\end{matrix} \right.} \right. & \left( {{II}{.4}} \right)\end{matrix}$

Regarding the nature of the problem, these harmonic functions are of thefollowing forms (e.g. for x-direction):

$\begin{matrix}{{h_{\Psi}\left( {\sqrt{k_{x}^{2} + \gamma_{x}^{2}}x} \right)} \sim \left\{ {{\begin{matrix}{\sin \left( {\sqrt{k_{x}^{2} + \gamma_{x}^{2}}x} \right)} \\{\cos \left( {\sqrt{k_{x}^{2} + \gamma_{x}^{2}}x} \right)} \\{e^{{\pm \sqrt{k_{x}^{2} + \gamma_{x}^{2}}}x}\mspace{59mu}}\end{matrix}{h_{~}\left( {\gamma_{x}x} \right)}} \sim \left\{ \begin{matrix}{\sin \left( {\gamma_{x}x} \right)} \\{\cos \left( {\gamma_{x}x} \right)} \\{e^{{\pm j}\; \gamma_{x}x}\mspace{25mu}}\end{matrix} \right.} \right.} & \left( {{II}{.5}} \right)\end{matrix}$

Substituting (II.5) in (II.1), one of the solutions for Equation (18)(referred to as the “elementary wave function”) is:

ψ_(n) =h _(Ψ)(√{square root over (k _(x) ²+γ_(x) ²)}x)h _(Ψ)(√{squareroot over (k _(y) ²+γ_(y) ²)}y)h _(Ψ)(√{square root over (k _(z) ²+γ_(z)²)}z)

φ_(n) =h _(Φ)(γ_(x) x)h _(Φ)(γ_(y) y)h _(Φ)(γ_(z) z)  (II.6)

However, as any two of these solutions are linearly independent (as per(II.3), among the three separation parameters (k_(x), k_(y), k_(z)) or(γ_(x), γ_(y), γ_(z)), only two of them are independent), the generalsolution to Equation (18) is the linear combination of (II.6) on the twoindependent functions (say (k_(x),k_(y)) or (γ_(x), γ_(y))) as

$\begin{matrix}{{\psi = {\sum\limits_{\sqrt{k_{x}^{2} + \gamma_{x}^{2}}}{\sum\limits_{\sqrt{k_{y}^{2} + \gamma_{y}^{2}}}{C_{1}\psi_{n}}}}}{\phi = {\sum\limits_{\gamma_{x}}{\sum\limits_{\gamma_{y}}{C_{2}\phi_{n}}}}}} & \left( {{II}{.6}} \right)\end{matrix}$

where C₁ and C₂ are two constants.

According to either Equation (10) or Equation (16), and depending on theboundary conditions of the problem, a specific form of themultiplication of φ and ψ in (II.5) gives the electric field or magneticflux density. Consequently, the argument of any type of thismultiplication must be (k_(x), k_(y), k_(z)). In other words, thefollowing general equation must be satisfied:

±γ_(x)±√{square root over (k _(x) ²+γ_(x) ²)}=±k _(x)

±γ_(y)±√{square root over (k _(y) ²+γ_(y) ²)}=±k _(y)

±γ_(z)±√{square root over (k _(z) ²+γ_(z) ²)}=±k _(z)  (II.3)

Solving the above equation for (γ_(x), γ_(y), γ_(z)) yields:

$\begin{matrix}{{{2k_{x}\gamma_{x}} = 0}{{2k_{y}\gamma_{y}} = 0}{{2k_{z}\gamma_{z}} = \left. 0\rightarrow\left\{ \begin{matrix}{{\left( {\gamma_{x},\gamma_{y},\gamma_{z}} \right) = 0}\mspace{256mu}} \\{{or}\mspace{400mu}} \\{{{\left( {k_{x},k_{y},k_{z}} \right) = 0},{{this}\mspace{14mu} {solution}\mspace{14mu} {violates}}}\mspace{20mu}} \\{{{the}\mspace{14mu} {wave}\mspace{14mu} {theory}\mspace{14mu} {which}\mspace{14mu} {requires}\mspace{14mu} k} \neq 0.}\end{matrix} \right. \right.}} & \left( {{II}{.4}} \right)\end{matrix}$

Thus, the only possible choice is (γ_(x), γ_(y), γ_(z))=0. This analysisalso reveals the time-varying nature of ψ and the static nature of φ.

Appendix III: Numerical Aspects in FEM

In order to solve the equations (18) via a FEM formulation, thecorresponding integral representation of these partial differentialequations (PDEs) are derived. This integral representation is called thefunctional of the PDE. As a general form, equation (28a) can bere-written as:

−∇² u+β ² u=0  (III.1)

where u is the scalar function i.e. v or ψ (as the two potentials to besolved) and β is −k².

The functional of equation (III.1) is derived as

$\begin{matrix}{{F(u)} = {{\frac{1}{2}\underset{V}{\int{\int\int}}\left( {{{\nabla u}.{\nabla u}} + {\beta \; u^{2}}} \right){dV}} + \underset{\underset{{contributes}\mspace{14mu} {only}\mspace{14mu} {on}\mspace{14mu} S}{}}{\frac{1}{2}\underset{S}{\int\int}\eta \mspace{14mu} u^{2}{dS}}}} & \left( {{III}{.2}} \right)\end{matrix}$

where V is the 3D volume of the problem domain and S is the surface onwhich the third kind boundary conditions are applied (e.g. thedielectric interfaces or the surrounding surface). According to thethird kind boundary conditions, η=h, q or w. The linear tetrahedralelement seen in FIG. 7(a) contains 4 nodes and is used in FEM as atypical element to discretize the problem domain. The value of theunknown function u is approximated at each linear tetrahedral element by

$\begin{matrix}{u_{m} = {\sum\limits_{i = 1}^{4}\; {N_{i}^{m}u_{i}^{m}}}} & \left( {{III}{.3}} \right)\end{matrix}$

Where N is called the interpolation function and u^(m) _(i) is the valueof the unknown function at ith node of the mth element. Theinterpolation function for linear tetrahedral is given as

$\begin{matrix}{N_{i}^{m} = {\frac{1}{6\Gamma}\left( {a_{i}^{m} + {b_{i}^{m}x} + {c_{i}^{m}y} + {d_{i}^{m}z}} \right)}} & \left( {{III}{.4}} \right)\end{matrix}$

for which the coefficients are

$\begin{matrix}{{{a_{i}^{m} = {\frac{1}{6\Gamma}{\begin{matrix}u_{1}^{m} & u_{2}^{m} & u_{3}^{m} & u_{4}^{m} \\x_{1}^{m} & x_{2}^{m} & x_{3}^{m} & x_{4}^{m} \\y_{1}^{m} & y_{2}^{m} & y_{3}^{m} & y_{4}^{m} \\z_{1}^{m} & z_{2}^{m} & z_{3}^{m} & z_{4}^{m}\end{matrix}}}};{b_{i}^{m} = {\frac{1}{6\Gamma}{\begin{matrix}1 & 1 & 1 & 1 \\u_{1}^{m} & u_{2}^{m} & u_{3}^{m} & u_{4}^{m} \\y_{1}^{m} & y_{2}^{m} & y_{3}^{m} & y_{4}^{m} \\z_{1}^{m} & z_{2}^{m} & z_{3}^{m} & z_{4}^{m}\end{matrix}}}}}{{c_{i}^{m} = {\frac{1}{6\Gamma}{\begin{matrix}1 & 1 & 1 & 1 \\x_{1}^{m} & x_{2}^{m} & x_{3}^{m} & x_{4}^{m} \\u_{1}^{m} & u_{2}^{m} & u_{3}^{m} & u_{4}^{m} \\z_{1}^{m} & z_{2}^{m} & z_{3}^{m} & z_{4}^{m}\end{matrix}}}};{d_{i}^{m} = {\frac{1}{6\Gamma}{\begin{matrix}1 & 1 & 1 & 1 \\x_{1}^{m} & x_{2}^{m} & x_{3}^{m} & x_{4}^{m} \\y_{1}^{m} & y_{2}^{m} & y_{3}^{m} & y_{4}^{m} \\u_{1}^{m} & u_{2}^{m} & u_{3}^{m} & u_{4}^{m}\end{matrix}}}}}{\Gamma = {\frac{1}{6}{\begin{matrix}1 & 1 & 1 & 1 \\x_{1}^{m} & x_{2}^{m} & x_{3}^{m} & x_{4}^{m} \\y_{1}^{m} & y_{2}^{m} & y_{3}^{m} & y_{4}^{m} \\z_{1}^{m} & z_{2}^{m} & z_{3}^{m} & z_{4}^{m}\end{matrix}}}}} & \left( {{III}{.5}} \right)\end{matrix}$

where (x_(i) ^(m), y_(i) ^(m), z_(i) ^(m)) are the coordinates of ithnode at mth element. Substituting (III.3) in (III.2) for each elementyields the functional at each element called F^(m). The Ritz method invariational principle says that the extremum point of F^(m) with respectto u^(m) _(i) is the solution of (III.1) at each element. Thus, to findthe extremum points, the derivate of F^(m) with respect to u^(m) _(i) isequated to zero at each element and the corresponding roots are found asthe solution. It reads

$\begin{matrix}{\frac{\partial{F\left( u_{i}^{m} \right)}}{\partial u_{i}^{m}} = {{{\sum\limits_{j = 1}^{4}\; {u_{j}^{m}\underset{V^{m}}{\int{\int\int}}\left( {{{\nabla N_{i}^{m}}.{\nabla N_{j}^{m}}} + {\beta \; N_{i}^{m}N_{j}^{m}}} \right){dV}}} + \underset{\underset{{contributes}\mspace{14mu} {only}\mspace{14mu} {on}\mspace{14mu} S}{}}{\sum\limits_{j = 1}^{3}\; {u_{j}^{m}\underset{S^{m}}{\int\int}\eta \mspace{14mu} N_{i}^{m}N_{j}^{m}{dS}}}} = 0}} & \left( {{III}{.6}} \right)\end{matrix}$

After writing equation (III.6) in a matrix form which gives the linearsystem of equations, we have

K_(4×4)U_(4×1)=O_(4×1)  (III.7)

where K and O are the coefficient and excitation matrices, respectively,and U is the matrix of unknowns. Their corresponding elements are

$\begin{matrix}{{K_{ij}^{m} = {{\underset{V^{m}}{\int{\int\int}}\left( {{{\nabla N_{i}^{m}}.{\nabla N_{j}^{m}}} + {\beta \; N_{i}^{m}N_{j}^{m}}} \right){dV}} + \underset{\underset{{contributes}\mspace{14mu} {only}\mspace{14mu} {on}\mspace{14mu} S}{}}{\underset{S^{m}}{\int\int}\eta \mspace{14mu} N_{i}^{m}N_{j}^{m}{dS}}}}{O_{ij}^{m} = {0\mspace{14mu} \left( {{except}\mspace{14mu} {the}\mspace{14mu} {Dirichlet}\mspace{14mu} {condition}} \right)}}} & \left( {{III}{.8}} \right)\end{matrix}$

The Galerkin's method is another technique to reach above results, butas it coincides with the Ritz method in equation (III.7) or (III.8), weexclude the discussion on this method.

As (III.3) is continuous, the continuity conditions in equations (23b)are satisfied automatically (u⁺=u⁻). Moreover, the continuity conditionson normal and tangential component is automatically satisfied splittingthe domain to subdomains where constitutive parameters (permittivity andpermeability) are continuous. Also, the surface integral in (III.2)includes the contribution of the third kind boundary conditions. Hence,the only remained condition to apply is the Dirichlet boundarycondition. To incorporate this condition on ith node, the followingchange in the corresponding matrix element of either K or O isnecessary.

K_(ii) ^(m)=a very large number, (say 10⁷⁰ as per [6])

O _(ii) ^(m) =f×(the above very large number)  (III.8)

When the system of equations is derived for each element, they areassembled (as seen in FIG. 7(b)) and the global system of equationsgiven bellow

U _(M×1)=(K _(M×M))⁻¹ O _(M×1)  (III.9)

is the solution of (III.1) where M is the total number of elements(3×10⁵ in our calculations). This is the general approach by which FEMsolves (III.1).

Appendix IV: Determining the Eigenvalue

In order to find the eigenvalue, the general forms of time-harmonicelectromagnetic waves are considered. The general travelling, one-sidedevanescent, or attenuated travelling solutions (according to thedifferent values of γ and √{square root over (2k²−γ²)}) for thepair-potentials described by (33) are:

$\begin{matrix}{{\psi = {{C_{0}e^{{- j}\sqrt{{2k^{2}} - \gamma^{2}}r}} + {C_{1}e^{{+ j}\sqrt{{2k^{2}} - \gamma^{2}}r}}}}{\phi = {{C_{2}e^{{- \gamma}\; r}} + {C_{3}e^{{+ \gamma}\; r}}}}} & \left( {{IV}{.1}} \right)\end{matrix}$

The general standing, two-sided evanescent, or localized standingsolutions (according to different values of γ and √{square root over(2k²−γ²)} for the pair-potentials described by (33) are:

ψ=C′ ₀ sin√{square root over (2k ²−γ²)}r)+C′ ₁ cos(√{square root over(2k ²−γ²)}r)

φ=C′ ₂ sinh(γr)+C′ ₃ cosh(γr)  (IV.2)

where C_(i) and C′_(i) are functions of space coordinates.

According to either (5) or (31), the multiplication of φ and ψ in either(IV.1) or (IV.2) gives the electric field or magnetic flux density.Consequently, this multiplication must have a wave number equal to jk.Hence, the following general equation must be satisfied

±γ±j√{square root over (2k ²−γ²)}=±jk

→±j√{square root over (2k ²−γ²)}=±jk∓γ  (IV.3)

where the left hand side of above equation is seen in the argument ofthe multiplication of φ and ψ in either (IV.1) or (IV.2). Solving aboveequation for γ yields

$\begin{matrix}{\left( {{\pm j}\sqrt{{2k^{2}} - \gamma^{2}}} \right)^{2} = {\left( {{\pm {jk}}\overset{\_}{+}\gamma} \right)^{2} = {\left. \left( {{jk} - \gamma} \right)^{2}\rightarrow{{{- 2}k^{2}} + \gamma^{2}} \right. = {\left. {{- k^{2}} + \gamma^{2} - {2{jk}\; \gamma}}\rightarrow{k^{2} - {2{jk}\; \gamma}} \right. = {\left. 0\rightarrow{k\left( {k - {2j\; \gamma}} \right)} \right. = \left. 0\rightarrow\left\{ \begin{matrix}{{\gamma = {{- \frac{j}{2}}k}}\mspace{315mu}} \\{{or}\mspace{391mu}} \\{{k = 0},{{this}\mspace{14mu} {solution}\mspace{14mu} {is}\mspace{14mu} {in}\mspace{14mu} {contradiction}}} \\{{{{with}\mspace{14mu} {wave}\mspace{14mu} {theory}\mspace{14mu} {requiring}\mspace{14mu} k} \neq 0.}\mspace{40mu}}\end{matrix} \right. \right.}}}}} & \left( {{IV}{.4}} \right)\end{matrix}$

Thus, the only possible choice is

$\gamma = {{- \frac{j}{2}}{k.}}$

1. A computational electromagnetics process for determiningelectromagnetic fields, the process being executed by at least oneprocessor of a data processing system, and including the steps of: (i)solving a pair of differential equations for only two variables, the twovariables representing a pair of scalar potentials; and (ii) determiningat least one of an electric field and a magnetic field from the pair ofscalar potentials determined at step (i).
 2. The computationalelectromagnetics process of claim 1, wherein step (ii) includes at leastone of: (i) determining the electric field {right arrow over (E)}according to:{right arrow over (E)}=jωψ∇φ; and (ii) determining the magnetic field{right arrow over (B)} according to:{right arrow over (B)}=∇φ×∇ψ; where φ and ψ represent the scalarpotentials, ω represents angular frequency, and j=√{square root over(−1)}.
 3. The computational electromagnetics process of claim 1, whereinthe differential equations are of the form:∇² ψ+k ²ψ=0∇²φ=0 where k=ω√{square root over (εμ)} is the wavenumber, ε representspermittivity, and μ represents permeability.
 4. The computationalelectromagnetics process of claim 1, wherein the differential equationsare of the form:∇² ψ+k ²ψ=0∇² v+k ² v=0 where k=ω√{square root over (εμ)} is the wavenumber, εrepresents permittivity, μ represents permeability, and v representsscalar electric potential, and wherein$v = {- {\frac{j\; {\omega\phi\psi}}{2}.}}$
 5. The computationalelectromagnetics process of claim 1, wherein the differential equationsare of the form: ∇²ψ + (2k² − γ²)ψ = 0 ∇²ϕ − γ²ϕ = 0 or∇²ψ + (2k² − γ²)ψ = 0${{\nabla^{2}v} + {k^{2}v}} = {- \frac{\rho + \rho_{j}}{ɛ}}$ where φand ψ represent the scalar potentials, ${\gamma = {{- \frac{j}{2}}k}},$k=ω√{square root over (εμ)} is the wavenumber, ε representspermittivity, μ represents permeability, v represents scalar electricpotential, and v=−jωφψ.
 6. A computational electromagnetics processexecuted by at least one processor of a data processing system, theprocess including at least one of: (i) determining an electric field{right arrow over (E)} according to:{right arrow over (E)}=jωψ∇φ; and (ii) determining a magnetic field{right arrow over (B)} according to:{right arrow over (B)}=∇φ×∇ψ; where ω represents angular frequency,j=√{square root over (−1)}, and φ and ψ represent scalar potentialssatisfying a pair of differential equations for only two variables, thetwo variables representing the scalar potentials.
 7. The computationalelectromagnetics process of claim 6, wherein the pair of differentialequations are of the form:∇² ψ+k ²ψ=0∇²φ=0 where k=ω√{square root over (εμ)} is the wavenumber, ε representspermittivity, and μ represents permeability.
 8. The computationalelectromagnetics process of claim 6, wherein the pair of differentialequations are of the form:∇² ψ+k ²ψ=0∇² v+k ² v=0 where k=ω√{square root over (εμ)} is the wavenumber, εrepresents permittivity, μ represents permeability, and v representsscalar electric potential, and wherein$v = {- {\frac{j\; {\omega\phi\psi}}{2}.}}$
 9. The computationalelectromagnetics process of claim 6, wherein the pair of differentialequations are of the form: ∇²ψ + (2k² − γ²)ψ = 0 ∇²ϕ − γ²ϕ = 0 or∇²ψ + (2k² − γ²)ψ = 0${{\nabla^{2}v} + {k^{2}v}} = {- \frac{\rho + \rho_{j}}{ɛ}}$ wherek=ω√{square root over (εμ)} is the wavenumber,${\gamma = {{- \frac{j}{2}}k}},$ ε represents permittivity, μrepresents permeability, v represents scalar electric potential, andv=−jωφψ.
 10. At least one computer-readable storage medium having storedthereon executable instructions that, when executed by at least oneprocessor of a data processing system, cause the at least one processorto execute the computational electromagnetics process of claim
 1. 11. Atleast one computer-readable storage medium having stored thereonexecutable instructions that, when executed by at least one processor ofa data processing system, cause the at least one processor to execute acomputational electromagnetics process, including the steps of: (i)solving a pair of differential equations for only two variables, the twovariables representing a pair of scalar potentials; and (ii) determiningat least one of an electric field and a magnetic field from the pair ofscalar potentials determined at step (i).
 12. The at least onecomputer-readable storage medium of claim 11, wherein step (ii) includesat least one of: (iii) determining the electric field {right arrow over(E)} according to:{right arrow over (E)}=jωψ∇φ; and (iv) determining the magnetic field{right arrow over (B)} according to:{right arrow over (B)}=∇φ×∇ψ; where φ and ψ represent the scalarpotentials, ω represents angular frequency, and j=√{square root over(−1)}.
 13. The at least one computer-readable storage medium of claim11, wherein the differential equations are of the form:∇² ψ+k ²ψ=0∇²φ=0 where k=ω√{square root over (εμ)} is the wavenumber, ε representspermittivity, and μ represents permeability.
 14. The at least onecomputer-readable storage medium of claim 11, wherein the differentialequations are of the form:∇² ψ+k ²ψ=0∇² v+k ² v=0 where k=ω√{square root over (εμ)} is the wavenumber, εrepresents permittivity, μ represents permeability, and v representsscalar electric potential, and wherein$v = {- {\frac{j\; {\omega\phi\psi}}{2}.}}$
 15. The at least onecomputer-readable storage medium of claim 11, wherein the differentialequations are of the form: ∇²ψ + (2k² − γ²)ψ = 0 ∇²ϕ − γ²ϕ = 0 or∇²ψ + (2k² − γ²)ψ = 0${{\nabla^{2}v} + {k^{2}v}} = {- \frac{\rho + \rho_{j}}{ɛ}}$ where φand ψ represent the scalar potentials, ${\gamma = {{- \frac{j}{2}}k}},$k=ω√{square root over (εμ)} is the wavenumber, ε representspermittivity, μ represents permeability, v represents scalar electricpotential, and v=−jωφψ.
 16. A computational electromagnetics systemhaving a memory and at least one processor configured to execute acomputational electromagnetics process, including the steps of: (i)solving a pair of differential equations for only two variables, the twovariables representing a pair of scalar potentials; and (ii) determiningat least one of an electric field and a magnetic field from the pair ofscalar potentials determined at step (i).
 17. The computationalelectromagnetics system of claim 16, wherein step (ii) includes at leastone of: (iii) determining the electric field {right arrow over (E)}according to:{right arrow over (E)}=jωψ∇φ; and (iv) determining the magnetic field{right arrow over (B)} according to:{right arrow over (B)}=∇φ×∇ψ; where φ and ψ represent the scalarpotentials, ω represents angular frequency, and j=√{square root over(−1)}.
 18. The computational electromagnetics system of claim 16,wherein the differential equations are of the form:∇² ψ+k ²ψ=0∇²φ=0 where k=ω√{square root over (εμ)} is the wavenumber, ε representspermittivity, and μ represents permeability.
 19. The computationalelectromagnetics system of claim 16, wherein the differential equationsare of the form:∇² ψ+k ²ψ=0∇² v+k ² v=0 where k=ω√{square root over (εμ)} is the wavenumber, εrepresents permittivity, μ represents permeability, and v representsscalar electric potential, and wherein$v = {- {\frac{j\; {\omega\phi\psi}}{2}.}}$
 20. The computationalelectromagnetics system of claim 16, wherein the differential equationsare of the form: ${\gamma = {{- \frac{j}{2}}k}},$ where φ and ψrepresent the scalar potentials, ∇²ψ + (2k² − γ²)ψ = 0∇²ϕ − γ²ϕ = 0 or ∇²ψ + (2k² − γ²)ψ = 0${{\nabla^{2}v} + {k^{2}v}} = {- \frac{\rho + \rho_{j}}{ɛ}}$k=ω√{square root over (εμ)} is the wavenumber, ε representspermittivity, μ represents permeability, v represents scalar electricpotential, and v=−jωφψ.